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ABSTRACT 


This  report  summarizes  the  results  of  a  study  on  nonlinear  phenomena 
in  plasms  with  the  min  emphasis  on  develoning  deterministic  theories  for 
plasm  turbulence  based  on  the  Ruelle-Takens  theory  of  strange  attractors 
and  on  the  interaction  of  solitons.  The  min  results  include  a  demonstra¬ 
tion  of  the  occurrence  of  turbulence-like  solutions  in  a  simplified  Zakharov 
model  for  Langmuir  turbulence  and  a  number  of  nonlinear  three -wave  inter¬ 
actions.  The  stability  of  mode -converted  lower-hybrid  solitons  is  also 
studied. 
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I.  INTRODUCTION 

Recently,  a  considerable  amount  of  work  has  been  performed  on  plasma 
turbulence.  This  was  partially  motivated  by  the  observance  of  various  forms 
of  turbulence  phenomena  in  plasma  heating  experiments  conducted  both  in  this 
country  and  abroad. 

Since  plasma  turbulence  is  primarily  due  to  the  nonlinear  excitation 
of  collective  oscillations  having  a  broad  spectrum  of  frequencies,  nonlinear 
plasma  models  must  be  used  in  the  development  of  theories  for  such  phenomena. 
The  classical  approach  to  the  theory  of  plasma  turbulence  is  to  regard  the 
phenomena  as  stochastic  processes.  Thus,  the  description  of  turbulence 
quantities  is  given  in  terms  of  their  statistical  averages  [1],[2].  In 
existing  theories  for  plasma  turbulence,  certain  basic  assumptions  are  in¬ 
troduced,  such  as  ergodicity  or  stationarity  of  the  processes.  The  validity 
of  these  assumptions  is  extremely  difficult,  if  not  impossible,  to  establish 
by  experiments.  Furthermore,  various  phenomena  at  the  onset  of  turbulence 
cannot  be  explained  by  statistical  theories. 

Motivated  by  such  deficiences  in  the  statistical  theories  of  turbu¬ 
lence,  attempts  have  been  made  to  develop  turbulence  theory  using  determin¬ 
istic  models.  In  1973,  Kingsep,  Rudakov  and  Sudan  [3]  suggested  that  strong 
Langmuir  turbulence  could  be  described  in  terms  of  a  system  of  interacting 
Langmuir  solitary  waves.  Subsequently,  a  considerable  amount  of  analytical 
and  numerical  work  has  been  done  based  on  this  idea  for  the  one -dimensional 
case  [4], [5].  It  was  suggested  that  a  similar  theory  could  be  developed  for 
the  multi -dimensional  case.  Goldman  and  Nicholson,  however,  showed  that 
multi -dimensional  spherically  symmetric  Langmuir  solitons,  under  adiabatic 
conditions  and  negligible  ion  inertia  effects,  are  unstable  [6].  Similar 
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conclusions  have  been  deduced  for  other  cases  [7] -[11].  Therefore,  further 
studies  are  necessary  before  definite  conclusions  can  be  drawn  from  the 
approach  suggested  by  Kingsep  et  al. 

In  early  1978,  under  the  support  of  AFOSR  Grant  No.  74-2662,  we  began 
to  explore  a  new  direction  for  developing  a  deterministic  theory  for  plasma 
turbulence.  This  approach  was  motivated  by  the  work  of  Lorenz  on  determin¬ 
istic  nonperiodic  flows  pertaining  to  atmospheric  turbulence  [1]  and  the 
theory  of  "strange  attractors"  proposed  by  Ruelle  and  Takens  [13] .  A  remark¬ 
able  feature  of  Lorenz's  work  is  his  demonstration  of  the  existence  of 
turbulence- like  solutions  of  a  single  nonlinear  third-order  ordinary  differ¬ 
ential  equation  with  quadratic  nonlinear  terms.  His  equation  is  a  highly 
simplified  model  of  an  incompressible  fluid  with  heat  transfer.  The  results 
suggest  the  possibility  of  modeling  atmospheric  turbulence  by  deterministic 
nonlinear  differential  equations. 

The  main  objective  of  this  study  is  to  develop  deterministic  theories 
for  plasma  turbulence  from  the  standpoint  of  the  Ruelle  and  Takens  theory 
of  strange  attractors  and  the  interaction  of  solitary  waves.  The  results 
of  this  study  will  be  summarized  categorically  in  the  subsequent  section. 
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II.  RESEARCH  SUMMARY 

In  the  study  of  plasma  turbulence  via  the  Ruelle-Takens  approach, 
attention  has  been  focused  on  Langmuir  turbulence  and  nonlinear  wave-wave 
interactions  in  plasmas.  The  study  of  plasma  turbulence  via  the  interaction 
of  solitary  waves  has  been  limited  to  mode -converted  lower-hybrid  solitons. 
The  results  are  sumnarized  below. 

A.  Nonperiodic  Oscillations  of  Langmuir  Naves 

In  1978,  under  the  support  of  AFOSR  Grant  No.  74-2662,  we  discovered 
that  single-mode  equations ,  derived  from  Zakharov's  model  for  Langmuir  turbu¬ 
lence  in  a  plasma  in  the  presence  of  an  external  spatially  homogeneous 
electric  field  oscillating  at  the  electron  plasma  frequency,  have  non¬ 
periodic  chaotic  solutions  whose  power  spectra  have  turbulence -like  features 
[14].  In  this  simplified  model,  phenomenological  damping  terms  are  intro¬ 
duced;  such  dissipative  terms  are  essential  for  the  existence  of  chaotic 
oscillations . 

During  the  period  covered  by  this  report,  we  explored  further  the 
observed  chaotic  oscillations  with  the  objective  of  achieving  some  under¬ 
standing  of  the  nature  of  the  onset  of  these  oscillations  and  their  quali¬ 
tative  behavior.  Our  first  results  consist  of  estimates  for  the  amplitudes 
of  the  chaotic  oscillations.  We  showed  that  the  magnitude  of  the  high- 
frequency  electric  field  oscillations,  about  some  real  scalar  multiple  of 
the  external  electric  field  E0,  satisfies  a  bound  which  grows  linearly  with 
||E0||  .  In  addition,  the  ion  density  oscillations  are  contained  in  a  ball 
whose  radius  grows  with  [| Eq || 2  .  These  bounds  are  consistent  with  the 
numerical  solutions.  The  results  were  combined  with  those  in  [14]  and 
published  recently  in  the  Journal  of  Mathematical  Physics  (see  Appendix  A) . 
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In  order  to  understand  the  onset  of  chaotic  oscillations,  we  determined 
numerically  the  Po incard  mappings  associated  with  various  hyperplanes  in 
the  system's  state  space.  It  was  found  that  the  mappings  are  quasi -one - 
dimensional,  but  their  structures  are  not  sufficiently  simple  to  permit  mean¬ 
ingful  analytical  studies.  Finally,  a  bifurcation  analysis  of  the  system  was 
performed  with  respect  to  the  damping  coefficients.  The  results  are  described 
in  the  report  [P9]  (see  Appendix  B). 


B.  Chaotic  Solutions  of  Resonant  Three-Wave  Interactions 
This  study  began  with  the  following  well-known  simplified  model  gov¬ 
erning  the  nonlinear  resonant  interaction  of  three  waves  with  both  linear 
growth  and  damping  terms: 


dAi/dt  =  YiA!  -  iVA2A2  , 
dA2/dt  =  -  y2A2  -  iVAjAj*  , 
dA3/dt  =  -  y3A3  -  iVAjA2*  , 


(1) 


where  the  Aj ' s  correspond  to  the  normalized  complex  wave  amplitudes ,  V 
is  a  real  coupling  coefficient  and  the  Yj ' s  are  positive  numbers.  This 
model  includes  a  number  of  important  resonant  wave -wave  interactions  in 
plasmas.  It  was  shown  that,  under  certain  conditions,  the  asymptotic 
behavior  of  the  solutions  of  (1)  for  t-*-°°  can  be  described  by  a  system 
of  three  real  ordinary  differential  equations.  Numerical  studies  showed 
that  the  reduced  system  exhibits  chaotic  oscillations  which  are  sensitive 
to  variations  in  the  initial  conditions.  Moreover,  the  autocovariance  func¬ 
tions  associated  with  the  chaotic  solutions  appear  to  decay  as  the  separa¬ 
tion  time  tends  to  00 .  This  suggests  that  a  turbulent  state  could  be 
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produced  through  the  nonlinear  interaction  of  a  linearly  unstable  wave  with 
two  linearly  damped  waves.  The  preliminary  results  of  this  study  are  de¬ 
scribed  in  the  report  [P7]  (see  Appendix  C)  published  in  November,  1979. 
Since  then,  further  studies  have  been  made  in  categorizing  general  three- 
wave  interactions  which  are  capable  of  producing  chaotic  oscillations.  The 
results  will  be  described  in  a  forthcoming  thesis  by  K.  Masui  [P8] . 

C.  Mode-Converted.  Lower-Hybrid  Soli  tons 

It  is  well  known  that  one -dimensional  Langmuir  solitons  are  modula- 
tionally  unstable  to  transverse  perturbations.  The  instability  causes  the 
solitons  to  collapse  and  form  bunches  resembling  a  turbulent  state.  Here, 
it  was  shown  that  the  planar  mode -converted  lower-hybrid  solitons  are  also 
unstable  to  transverse  perturbations.  Moreover,  two  distinct  classes  of 
modes  are  found.  In  the  lowest  order,  the  stability  of  each  class  is  re¬ 
lated  to  a  particular  term  of  the  equation  governing  the  nonlinear  evolution 
of  these  waves.  The  details  are  described  in  the  report  [P5]  (see  Appen¬ 
dix  D). 


D.  Chaotic  Cscillations  in  Bilinear  Systems 

It  was  observed  that  the  mathematical  models  of  may  physical  processes, 
including  certain  nonlinear  wave -wave  interactions  and  the  Lorenz  equation, 
can  be  written  in  the  form:  dx/dt  *  Ax  +  f(x)  ,  where  A  is  a  real  n><n 
matrix  and  f  is  a  bilinear  function  of  x  .  During  1979,  chaotic  oscilla¬ 
tions  in  this  class  of  systems  were  explored  from  the  viewpoint  of  feedback 
control  theory.  A  simple  and  sufficient  condition  for  the  nonexistence  of 
such  oscillations  was  derived.  Ellipsoidal  bounds  for  the  amplitudes  of 
chaotic  oscillations  were  also  obtained.  The  results  are  applicable  to  the 
Lorenz  equation  and  other  equations  which  are  known  to  have  chaotic  oscilla¬ 
tions.  The  details  are  given  in  report  [P6]  (see  Appendix  E) . 
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III.  CONCLUDING  REMARKS 

Since  the  initiation  of  this  study,  there  appears  to  be  a  growing 
interest  in  the  theory  of  chaotic  oscillations  in  deterministic  nonlinear 
dynamical  systems  and  its  application  to  turbulence,  as  proposed  by  Ruelle 
and  Takens .  During  the  period  covered  by  this  report,  the  following  confer¬ 
ences  and  symposia  were  held  both  in  the  US  and  abroad  at  which  papers  on 
chaotic  oscillations  were  presented: 

(1)  International  Workshop  on  Intrinsic  Stochasticity  in  Plasmas, 

June  17-23,  1973,  Institut  d'Etudes  Scientifiques  de  Car gese, 
Corse,  France. 

(2)  Annual  Meeting  of  the  Plasma  Physics  Division,  American  Physical 
Society,  November  12-16,  1979,  Boston. 

(3)  Engineering  Foundation  Converence  on  New  Approaches  to  Nonlinear 
Problems  in  Dynamics,  December  9-14,  1979,  Asilomar,  California. 

(4)  International  Conference  on  Nonlinear  Dynamics,  December  17-21, 
1979,  New  York  City. 

So  far,  chaotic  oscillations  in  nonlinear  systems  have  been  discovered 
mostly  through  bifurcation  analysis  and  numerical  experimentation.  There 
are  no  general  mathematical  conditions  of  ensuring  the  existence  of  chaotic 
oscillations  in  systems  described  by  ordinary  or  partial  differential  equa¬ 
tions.  Moreover,  it  is  very  difficult  to  verify  that  a  numerically  computed 
apparently  chaotic  solution  is  indeed  chaotic  in  some  sense.  Until  these 
basic  questions  are  resolved,  there  remains  a  gap  between  the  mathematical 
theory  of  chaotic  oscillations  and  turbulence  phenomena  as  observed  in  the 
physical  world. 
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Nonperiodic  oscillations  of  Langmuir  waves 

P.  K.  C.  Wang 

School  of  Engineering  and  Applied  Science.  University  of  California.  Los  Angeles.  California  90024 

(Received  21  December  1978;  accepted  for  publication  15  October  1979) 

It  is  shown  that  the  single-mode  equations  derived  from  Zakharov’s  model  for  Langmuir 
turbulence  in  a  plasma  in  the  presence  of  an  external  spatially  homogeneous  electric  field 
oscillating  at  the  electron  plasma  frequency  has  nonperiodic  chaotic  solutions  whose  power 
spectra  have  turbulence-like  features.  Bounds  for  these  chaotic  solutions  are  derived.  Typical 
numerical  results  are  presented  for  the  one-dimensional  case. 


1.  INTRODUCTION 


It  has  been  observed  that  certain  nonlinear  ordinary 
differential  equations  have  chaotic  or  turbulence-like  solu¬ 
tions.1'1  A  simple  example  is  the  Lorenz  model  for  thermal 
convection  in  a  fluid  layer.1  Recently,  Ruelle  and  Takens 
proposed  that  fluid  turbulence  can  be  mathematically  char¬ 
acterized  by  this  class  of  solutions  whose  trajectories  in  the 
state  space  are  attracted  to  a  nonempty  set  ("strange  attrac¬ 
tor’’)  which  is  neither  an  equilibrium  set  nor  a  periodic  or¬ 
bit.’  On  this  set,  the  trajectories  exhibit  chaotic  oscillations. 
Moreover,  they  are  sensitive  to  variations  in  the  initial  condi¬ 
tions.  Here,  we  shall  demonstrate  that  the  single-mode  equa¬ 
tions  derived  from  the  Zakharov’s  model  for  Langmuir  tur¬ 
bulence  in  a  plasma  have  nonperiodic  chaotic  solutions. 

We  begin  with  the  following  dimensionless  form  of 
Zakharov’s  equations  describing  the  nonlinear  interaction  of 
high-frequency  electron  oscillations  with  an  ion  fluid  in  the 
presence  of  an  external  spatially  homogeneous  electric  field 
oscillating  at  the  electron  plasma  frequency  to  p  :M 


V-[i  — -  +V2E-«(E  +  ^,) 


(1) 


V2”  =  V2[|E|J  +  E(,-(E  7-  E*)] ,  (2) 

<9/ 2 


where  i  =*  V  —  1;  (•)*  denotes  complex  conjugation,  and  a  b 
the  usual  scalar  product  of  two  real  or  complex  vectors  a  and 

b.  E  =  {£,  . £v)  is  the  complex  amplitude  of  the  high- 

frequency  electric  field  $  given  by 

&(t,  x)=«Re[E(r,x)exp(-is»J,r)];  (3) 

and  n  is  a  real  quantity  corresponding  to  the  low-frequency 
perturbation  in  the  ion  density  from  its  constant  equilibrium 
value  n0  .  The  units  of  time  r,  spatial  coordinates  x  =  (x, , 
...,xv),  electric  fields  Eo  and  E,  and  ion  density  perturbation 
are  respectively,  3/(2 ouu,),  (3/2 XD  ,  [(64/ 3)77710111, 
c2  ],/* ,  and  (4/3 )an0  ,  where  a  is  the  electron-ion  mass  ratio 
m./m, ,  Aa  the  Debye  length,  and  c,  is  the  ion  acoustic 
speed.  Here,  Eo  *  (£0 . ,£0/v)  is  a  real  constant  vector  cor¬ 

responding  to  the  normalized  amplitude  of  the  external  elec¬ 
tric  field.  It  is  of  interest  to  determine  the  behavior  of  the 
solutions  of  (1)  and  (2)  (with  appropriate  damping  terms)  as 
a  function  of  the  parameter  Eo  ,  in  particular,  the  existence 
of  turbulence-like  solutions  for  some  Eo  ,  and  the  onset  of 
such  solutions  as  Eq  tends  to  some  threshold  values. 


2.  SIMPLIFIED  MODEL 

Let  the  spatial  domain  fl  be  a  bounded  open  subset  of 
the  ,V-dimensionai  Euclidean  space  R'v,  and  L  2(/? )  denote 
the  Hilbert  space  of  real  square-integrable  functions  defined 
on  fl  with  inner  product  ( u.u )  =  5 „  u(x)v(x)  d(l.  Let  {dj. } 
be  a  countable  orthonormal  basis  for  L  \fl ).  We  seek  solu¬ 
tions  to  (1)  and  (2)  in  the  form: 


E(r,x)  =  £  Et  (r  V „  (x),  n(t,x)  *  „  (t  )4k  (x) .  (4) 

If  the  boundary  of  12  is  sufficiently  smooth,  then  the 
Lapiacian  with  suitable  homogeneous  boundary  conditions 
is  a  negative  operator  with  a  countable  point  spectrum.  We 
may  take  <bk  to  be  the  orthonormalized  eigenfunction  of  V2 
corresponding  to  the  eigenvalue  Ak  —  —  ft\  In  this  case, 
we  may  substitute  (4)  into  (1)  and  (2),  multiply  both  sides  of 
the  equations  by  i„  (x) ,  and  integrate  over  12  to  give  a  coun¬ 
tably  infinite  system  of  ordinary  differential  equations  for 
E„  and  nm : 

i — — - =  fl^Eo  4-  a„kk  rttEi,  ,  (5) 

at  k  *• 


£V 

dt1 


■+n'mn„ 


=  — /»*  Eo-fE^  +E;)+  £££,**  E4-Et.  ,  (6) 

where  E^E*.  =  'ZjEkjE\-1  and 

aw  =  I*  d*(x*Mx>Mx)<//2.  (7) 

Jn 


(8) 


0mkk  =  f  V2  [tk  (*>**  •  (*)  ]<t>m  (*)  df2 . 

Jn 

By  retaining  only  the  terms  involving  E„  and  nm  in  (5) 
and  (6),  we  obtain  the  following  simplified  equations  for  a 
single  mode  m : 
d  E 

i— --(Ml,  -ty.)Em  =nm(E„+awE^),  (9) 

ut 

dln” dnm±u.-n 


=  -qi.Eo-fE^  +E*W)+  Em  |E„.  |2 , 


(10) 


where 
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am  —  f  <bl(x)dn,  pm  =  f  v2[<**:(*)R.(*)d/2. 

Jn  Jn 

(ID 

Also,  we  have  added  the  phenomenological  damping  coeffi¬ 
cients  ym  and  r„  .  They  may  represent  Landau  damping  of 
the  high-  and  low-frequency  waves.  In  what  follows,  we  shall 
analyze  the  behavior  of  the  solutions  of  (9)  and  (10)  as  Eg 
varies.  For  brevity,  the  subscript  m  in  (9)  and  (10)  will  be 
omitted  in  places  where  ambiguity  does  not  arise. 

A  starting  point  for  searching  the  strange  attractor  or 
chaotic  solutions  is  to  study  the  nature  of  the  equilibrium 
points  as  Eg  varies.  It  is  known  that  chaotic  solutions  could 
arise  after  finite  number  of  Hopf  bifurcations,’  therefore  we 
shall  establish  the  existence  of  Hopf  bifurcation  points. 

3.  EQUILIBRIUM  POINTS 

Consider  the  following  equations  for  determining  the 
equilibrium  points  of  (9)  and  (10)  for  any  given  Eg : 

(iy  -  n2)E  =  /i(Eg  +  aE)  ,  (12) 

H2n  =  — /rJEg-(E  -f  E*)  4-  0  |E|2  .  (13) 

Using  (1 3)  to  eliminate  n  in  (12),  we  obtain  an  equation 
for  E: 

(iy  —  /i2)E  =  [  —  Eg-(E  +  E*)  +  0n  ~2  |E|*](Eg  +aE). 

(14) 

Evidently,  since  Eg  is  a  real  TV-dimensional  vector  and  n 
is  real,  a  solution  of  ( 14)  must  be  a  complex  scalar  multiple  of 
Eg  (i.e.,  E  —  £Eg  for  some  I-  =  £K  +  i£,).  Thus,  the  solution 
of  (14)  reduces  to  finding  £.  Note  that  E  =  0  is  a  solution  of 
(14)  for  any  Eg  .  Substituting  E  =  £  Eg  into  (14)  leads  to  the 
following  equations  for  £„  and  £, : 

-  (m2£k  +  y£i) 

=  +  £  /)  —  2£a  )(1  +  &£k  )||Eo  l|J,  (15) 

Hr  -  - *(i  i  +  £  ?)  -  2ft,  )ag,  ||3,  II2, 

(16) 

and 

«  =  ifif* -  *  +  £  ?)  -  Ur  }||Eo  r,  (17) 

where  ||Eo||2  =  Eg-Eg  . 

Dividing  (15)  by  (16)  gives 

£1 i+£J  =  («r)"V£/- tf*).  (18) 

which  implies  that  a  solution  must  lie  on  the  circle: 

Now,  we  substitute  (18)  into  (16)  and  solve  for  £R  in  I 


terms  of  g , : 

Sm  =S'(m'-+  fir-'  IIEgiri/) 

X  [y  +  (2a  +  0p~  :)||Eg  ||2>/  ]  ~  '  ■  (20) 

Finally,  using  (20)  to  eliminate  g  „  in  ( 1 8)  leads  to  the  follow¬ 
ing  quadratic  equation  for  g, : 

A£2,+B£,+C  =  0,  (21) 

where 

A  =  [(0/y)2  +  <52]||Eg||\  (22) 

B  =  2(ji20y~'  +  y5)||Eg  ||2  -  V<5y  -  '||Eg  ||4 .  (23) 

C  =  fi*  +  y2  —  2^t2||Eg  ||2  (24) 

5  =  (2a  +  0fi~2) .  (25) 

If  B 2  -  AC>0  or 

||Eg ir  -4 a  08-  2||Eg  II2  -  (lay/5)2>0  ,  (26) 

then  (21)  has  real  roots  given  explicitly  by 
fl=  srn2  fur.  ip _  (M20  +  fS) 

(02 + sy)^ ii2  r*!l  *i2 

±[||^ir-4a^-2||Eg||2-(^)2],/J)  (27) 

For  ||Eo  ||2>0,  condition  (26)  is  satisfied  if  and  only  if 
||Eg||2>E^. ^2a<5-2[  0  +  (02  +  ^y2)1'2] .  (28) 

Thus,  we  conclude  that  for  0<  ||Eg  ||  <  Ex  ,  the  origin 
(E,n,h)  =  (0  +  jO,  0,0)  is  the  only  equilibrium  state  of  system 
(9)  and  (10),  where  n  denotes  dn/dt.  When  ||Eg  ||  =  Ex  ,  a 
new  equilibrium  state  (E,/i,h)  =  (Re(E)  +  iIm(E),n,0) 
emerges,  where Re(E)  =  £*Eg,Im(E)  =  £*Eo  with £,  given 
by 

6  _  tyn1  (Fi  _  tffi+r1®)  (29) 

(02  +  s2y2)E2oc  \  Sfi2  )’ 

and  S*  <n  arc  given  by  (20)  and  (17)  respectively.  As  ||Eg  || 
increases  from  Ex  ,  the  foregoing  nonzero  equilibrium  state 
bifurcates  into  two  distinct  equilibrium  states  (£  £  E, 

+  i£  ,*E<,  ,n  *  ,0)  and  (£  *~  Eg  +  j£,-Eg  ,n  “  ,0)  where  £ 
are  given  by  (27),  whose  corresponding  £  and  n±  are  de¬ 
termined  respectively  by  (20)  and  (17).  We  note  that  £K  and 
£,  depend  on  ||Eg  ||2.  Also,  the  coefficient  C  defined  by  (24) 
vanishes  when 

||Eg||2  =  £^  =(//*  +  f)/{2fi2) .  (30) 

Consequently,  £ also  vanishes.  Thus,  in  this  case,  the  equi¬ 
librium  set  consists  of  the  origin  and  the  point  [(£  « 

+  i£  *  )Eg  ,n  +  ,0] .  When  ||Eg  ||  increases  from  Ex  ,  we 
have  again  three  distinct  equilibrium  points. 


4.  STABILITY  OF  EQUILIBRIUM 

Let  Eg  «*  Re(E),  E,  *  Im(E),  and  z  denote  the  2(TV  +  I)-dimensional  real  vector  (nrn,E„,E1)r,  where  (-)r denotes 
transposition.  We  rewrite  (9)  and  (10)  in  the  form: 


at 


-fn  -2Th-lti2  Eg-E*  +  0  (||E„||2  +  \\E,\\2) 
-  yE*  +(+i2  +  an)E, 

,  -  yE,  -  (fi2  +  an)E„  -  nEg 
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Letz,.  =  (fl‘,0,Ei  ,E))be  an  equilibrium  point  of  (3  l)as  given  in  Sec.  3,  and<5z(r  )^z(i )  —  zr  .  We  consider  the  following 
linearized  system  of  (3 1)  about  z. : 

~-J/( z.;Eo)6z.  (32) 

dt 

where  i  f  (zf;Eo)  is  the  Jacobian  matrix  of  f  at  z,  given  by 


J  /(zc;Ea)  — 


1 


os 


—  A* 
a£^ 

■(Eo  -t-aE4*) 


-2T  2(/?E*  — /<2Eo)r 

o  -Kv 

0  -  (^2  +  an^Iy 


oc 

20(E))r 
Ou2  +  anOIv 
-  rlv 


(33) 


where  0,v  and  Iv  are  the  iV-dimensional  zero  vector  and  N  X.N  identity  matrix  respectively.  It  can  be  shown  (see  Appendix) 
that  the  characteristic  polynomial  of  J/(z,;E0)  is  given  by 

=  (O'  +  'O2  +  (M*  +  ane)i]N(Az  +2TA.  +  ja2)  +2  [(y  +  A  )2  +  (fi2  +  an02l'v_l 

xfor2  +  «/»*){( 0E'k  -^2EoHEo+Ei)  +  a)3||E5!!2}+(y  +  /lKaM2+  (3)Eo-E',].  (34) 

For  the  case  where  z,  =*  0,  the  above  expression  reduces  to 

det^flhEJ-A  I*,*,,  ]  =  l(r  +  'D2  +mT~]  [(/l 2  4-2 H.  +  /){(y  +  A  )2  +  M*}  -2M4!|Eo||2]  (35) 

Evidently,  when  Eq  =*  0,  the  spectrum  of  J/  (0;0)  is  given  by  {  —  Y  ±  i/i2  (multiplicity  N),  —  T  ±  (r 1  —  ji2)u2]  .which 
implies  the  asymptotic  stability  of  the  origin  for  yJT  >  0.  For  Eq  ^  0,  the  eigen  values  A  *  -  y  ±  ifi1  [multiplicity  (/V  - 1 )] 
remain  invariant,  while  the  remaining  eigenvalues  are  roots  of  the  quartic  equation: 

A 4  +2(y  +  r)A3  +(/z*+fii  +4y/>  r)A 2  +2  [r(f  +fi*)  +  y/i2]A  -^[(y2  +  /r‘)  — 2 /a2||Eo (|2]  =  0 .  (36) 

Obviously,  J {  (OjEq)  has  a  zero  eigenvalue  when  ||E(,  |j2  =  E2X  +  y2)/^2)  .  This  coincides  with  condition  (30)  for 
which  one  of  the  equilibrium  points  returns  to  the  origin.  It  can  be  readily  shown  by  using  Routh’s  criterion’  that  the  origin 
becomes  unstable  when  ||E<,  ||  >  Ex  .  In  fact,  (36)  has  only  one  unstable  root  and  it  is  real  and  positive.  Thus,  the  origin  has  a 
saddle  point  structure  in  a  two-dimensional  manifold.  So  we  conclude  that  Hopf  bifurcation  cannot  occur  at  the  origin  for  any 
value  of  ||  Eo  || . 

For  the  case  where  ||E„  ,  there  exist  nonzero  equilibrium  states  z,  which  depend  on  ||Eo  || .  We  observe  from  (34) 

that  A  =*  —  y  ±  i </t2  +  an *)  are  stable  eigenvalues  of  J,  (z,  ;Eo)  with  a  multiplicity  of  (JV  — 1).  Tie  remaining  eigenvalues  are 
given  by  the  roots  of  the  quartic  equation: 

A  4  +  a3  A 3  +  u,  (HEolDA  2  +  a,  (HE,,  ||)A  +  a0(  HE,  II)  =  0 .  (37) 

where 


a3  =  2 (y  +  O,  fljdlE,, ||)  »  y*  +  n2  +4 rr  +  (fi2  +  an*)2 , 

(IIEJ)  =  2{r  [y2  +  (m2  +  a«*)2]  +  rfi2+(aj*2  +  0)Eo-E)},  (38) 

a0(||Eo||)  =M2ty2  +  ^2  +  <m')2]  +  2{<>*2  +  E'„  EoKEo  +  aE'*) 

+  a0||E)||2]  +  y(a#t2  +  ^)E0-E'/}, 

where  Ei  ,EJ ,  and  n'  (given  in  Sec.  3)  depend  on  j|E<,  ||  .To  determine  the  value  of  ||Eo||  for  which  Hopfbifurcationoccurs.it 
is  necessary  to  determine  the  existence  of  purely  imaginary  roots  of  (37)  for  some  value  of  j|Eo  || .  From  Routh’s  criterion,  we 
can  deduce  that  if 

a,  ai  (||Eo  ||)  >  a.  (||Eo  II)  (39) 


a,(|l^||)[a3fl1(||^||)-ai(||E0||)]  =a2a0(||E0||),  (40) 

then  (38)  has  a  pair  of  purely  imaginary  roots  given  by  A  =  ±  [ajajflIEo  ||)  —  a,  (j|Eo  ||)]},/2i . 

Let  Eoh  be  the  value  of  ||E,  ||  such  that  both  (39)  and  (40)  are  satisfied;  and  A'  ±  (||Eo  ||)  =*  A*(j|Eo  ||)  ±  iX, (|jEo  ||)  be  the 
roots  of  (37)  such  that  A  „  (E0H)  =  0.  By  a  lengthy  but  straightforward  computation,  it  can  be  shown  that/  R  ,  the  derivative  of 
X„  with  respect  to  the  parameter  ||Eo|| ,  is  given  by 

x  ■*  (||Eo  II)  =  WflodlEo  II)  +2  a,  (HE,  ||V,  (flE,  ||)  -  a,a2  (||Eo  ||)a,  (||Eo  ||) 

-  a3a2(||E0  fly,  (||Eo  ||)]{2a3  [a,«,  (||Eo  ||)  +  a*  ([|Eo  ||)  -4  a0(||E^  ||)]}  -•  ,  (41) 

where  a)  denotes  the  derivative  of  a,  with  respect  to  ||E^|| .  For  Hopf  bifurcation, 10  A  „  (£ow)  >  0 .  Due  to  the  complicated 
dependence  of  a2  .a, ,  and  a0  on  ([E^  [( ,  it  is  difficult  to  determine  the  threshold  values  of  ||Eo  ||  for  Hopf  bifurcation.  We  shall 
resort  to  numerical  computation  at  this  point. 
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FIG.  1.  Locus  of  equilibrium  electric  8eld£'>c£«  +iE',  with  El  u  a 
parameter.  (£,*  «•  $  *  £^  and  £  "  =*  £  ~  are  denoted  by  solid  dots  and 
circles  respectively. 


We  note  here  that  the  Hopf  bifurcation  problem  for  (9) 
and  (10)  with  spatial  dimension  TV  or  dim(E)>  1  can  be  com¬ 
pletely  studied  by  considering  only  (37)  which  is  the  charac¬ 
teristic  equation  for  the  case  with  N  —  1.  Since  for  />  0  and 
r>  0,  the  additional  eigenvalues^  =  —  A.  ±  i(p.2  +  an1)  for 
Af>  1  are  stable,  and  a2  ,fl,  (||Eo  H)  >  0  for  all  ||Eo  ||,  the  di¬ 
mension  of  the  unstable  manifold  associated  with  a  nonzero 
equilibrium  state  is  at  most  three. 

5.  BOUNDS  FOR  CHAOTIC  OSCILLATIONS 

The  existence  of  chaotic  oscillations  depends  on  the 
manner  in  which  the  stable  and  unstable  manifolds  associat¬ 
ed  with  the  equilibrium  points  intersect  with  each  other.  At 
present,  there  are  no  readily  verifiable  analytical  sufficient 
conditions  for  the  existence  of  chaotic  solutions  for  finite 
dimensional  systems  of  ordinary  differential  equations. 
Here,  we  assume  the  existence  of  chaotic  oscillations  and 
proceed  to  derive  bounds  for  their  amplitudes,  thus  provid¬ 
ing  estimates  for  the  size  of  the  invariant  manifold  generated 
by  the  chaotic  oscillations. 

Fust,  we  shall  make  use  of  a  function  V  of  the  form: 


FIG.  2.  Locus  of  equilibrium  ion  densititiei  n  *  and  n  *  with  £  a  is  t 
parameter. 
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FIG.  3.  Locus  of  ihe  eigenvalues  of  1,  ((>,£„)  with  £■  as  a  parameter. 

V(E)=  |E-E|V2,  (42) 

to  estimate  the  magnitude  of  chaotic  oscillations  of  the  elec¬ 
tric  field  where  E  =  ReE  +  /ImE  is  to  be  determined.  By 
direct  computation: 

-rGlE*  +  (i)ly-'(M1  +  a«)IinE-ReE]||J 

+  ||E,  +  (J)  [Y <M2  +  an)  ReE  —  ImE  +  y_lnE0  ]||2 
-  Ullr-'  (M2  +  an)ImE  -  ReEf 
+  Hr <M2  +  crn)ReE  —  ImE  +  y *'  *E„  ||2  ]> .  (43) 

If  we  set  ReE  =  —  a  _l  Eq  and  ImE  =  0,  then  (43)  reduces 
to 

~T~—  ~  rOlE*  +  (2a) -1  Eo ||J  +  ||E,  -m%/(2 ar)||J 

at 


FIG.  4.  Locus  of  the  eigenvalues  of  J,  (s'  :£,)  with  £„  as  a  parameter. 
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FIG.  5  Locus  of  the  eigenvalues  of  J ,  (x  ;£,)  wuh  £,;  as  a  parameter. 

-(2a)-2(l+/iV-2)||E0||2} 

=  —  f[  |E  —  E,  |2  —  (2a)  _2(1  +/i4j/-2)||E0||2] ,  (44) 

where  E,  =  —  (2 a)  Eo  +  //a2E0/(2a><).  We  note  that 
with  the  foregoing  choice  of  E,  n  does  not  appear  in  (44). 
Moreover,  for  y  >  0,  d  V /dt  <  0  at  any  point  E  exterior  to  the 

set  =■  =  {E:  j  E  -  E,  |  <(2a)  - 1  (1  +  mV  >‘/2||Eo  ||> .  Let 

■=«=*( E:  V(E) » J|£  +  Since  for  any  <5>0, 

■=*  is  a  bail  in  E-space  centered  about  the  point  —  a  Eq 
with  radius  V2S,  it  is  possible  to  select  a  <5  such  that  iCl{  . 
In  fact,  elementary  geometric  considerations  show  that  the 
smallest  <5  having  the  foregoing  inclusion  property  is  given 
by<5  =  (i /2a)-'  (1  +  MV'2)l/2||Ea|| .  Evidently, 
dV/dt<  0  at  any  point  E  exterior  toSs  .  This  implies  that  for 
a  solution  of  (31)  initiated  from  any  point 
z(0)  =  (n(0),/i(0),E*  (0),E,  (0))  at  t  =  0  with  E(0)  =  E*(0) 
+  iE,  (0)  exterior  to  E&  ,  its  corresponding  E(f ),  t  >  0  either 
eventually  enters  jsg  at  some  finite  time  t ,  >  0  and  remains  in 
£3  for  all  r  >f,  or  tends  to  £3  as  t  -*  00 .  Clearly,  £3  con- 


25.  cc  ?^.  ;o 


' 
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FIG.  6.  £(;>|  and  *(/ )  vs  lime  t  corresponding  to  the  solution  of  ( 53)  and 
(54)  with  m  I  .£»  r/vW.  T  -  ?.0.  y,  -  l.0,£-  -  1.625;  initial  data: 
£<0)  -  1.426  *  0.5071/.  n(0)  -  -  40.75  and  n(0)  »  34.58. 
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-6.00  -3.00  0.00  3.00  5.30 
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FIG.  7  Projection  of  the  trajectory  of  (53)  and  (54)  (with  parameters  as 
given  in  Fig.  6)  onto  the  (£,  ,£,)-plane. 


tains  all  the  points  E(f )  along  any  chaotic  solution  of  (3 1) 
when  it  exists. 

Next,  we  derive  a  bound  for  the  magnitude  of  ion  densi¬ 
ty  oscillations.  Let  h  =  n  +  fi20  ||EJ2.  We  can  rewrite 
the  first  two  equations  in  (31)  as 

/L  V 

^-na^r+tO^lE-M2^-^!2]7-,  (45) 

where -/F"  =  (n,/i)rand 


Given  *1^(0),  the  initial  data  for-F"  at  t  =  0,  (45)  is  equiv¬ 
alent  to  the  integral  equation: 

^V{t )  =  [expaaff  ]^F*(0)  +  J*  exp^f  -  r)] 

x[0,^|E(r)-M2^-'Eo|J]rdr.  (47) 

Thus, 

|pr(f)||<||exp^f  |||pn0)||  +  J'  |  0  |  !|expj/ (t  -  r)|| 

x|E(r)-M2/ff  _,Eo|2rfr.  °  (48) 

For  m2>  r  >  0,  we  can  find  a  constant  #  >  0  such  that 
||exp  (jafr )  j|<#  exp(  —  A ).  Also,  we  have  already  estab¬ 
lished  that  along  any  chaotic  solution  E(r  )C.=3  or  |  E(f ) 

+  a-lEo  |<<5  =  (\/2a)-l(l  +mV“2)i/2||Eo||  for  all  t. 
Hence, 

|E(r)  — m2£  Eq  | 

=  |E(r)  +  a  _l  Eq  —  (a-1  +  fi20~')Eo\ 

<|E(r)  +  a-'E0|  +  |a-'  +  ^ 0 |||E^|| 

<S+  |a"'  +m2(S~'  |||EoH<^||Eo||,  (49) 

where  V  =  (3/ 2a) (1  +  mV~2)i/2  4-  |a-'  +  m2 0  |. 

It  follows  from  (48)  that 

l|.r(oi|<^{!F|j?ii|Eo||2r-'  +[||*«»|| 

-^l^|||Eo||2r-']exp(-n)} 


402 


J.  Matt).  Ptiy»„  Vol.  21 ,  No.  2,  February  1 980 


P.K.C.  Wang 


402 


FIG.  8.  Solution  of  (S3)  and  (54)  with  El  -  2.649, 

£(0)=  -3.07656  +  4.22328/. n(0)«  -8.38056,6(0)  =  -8.25344,  and 
other  parameters  as  given  in  Fig.  6. 


<#max{||z(0)||,!P  |  0  |  ||E0||2/'  -1 } ,  (50) 

for  aU  t  >0  and  E(0)e=V  When  ||z(0)||<^|  0  |  !|E0||2r-1, 
we  have  0  |  ||E0|j2/''1  for  all  t  >0,  a  bound 

which  is  independent  of  z(0). 


6.  ONE-DIMENSIONAL  EXAMPLE 


Consider  the  case  where  N  *  1  with  a  bounded  spatial 
domain  l J  =  (OX  ].  Assuming  that  both  E and  n  vanish  at  the 
boundary  points  x  =  0  and  x  *  L,  we  can  take  6t  (x) 

*  (2 /L  )[,1i\n{kirx/L ),  fik  ==  kir/L,  k  =  1,2—.  For  this 
case,  the  coefficients  am  and  0m  defined  in  (1 1)  become 


r 


{2/L)in%m\mirx/L)dx 


0  for  m  even, 

(2VL )l/2/(3m5r)  for  m  odd,’ 


(51) 


0m  =  f  (  —  )  [im:(rmrx/L  )]sin(mirx/L  )  dx 

Jo  \  L  )  dxr 

f  0  for  m  even. 

I  —  8mir(2/£  5) 1 :/3  for  m  odd. 

Thus,  for  an  odd  integer  m,  Eqs.  (9)  and  (10)  have  the 
form: 


vi. 


o 

a 

d. 


74 - 1 - 1 - 1 - 1 - 1 

-15.00  -10.00  -S.00  0.00  5.00  10.00 

_ _  REtE) 

FIG.  9.  Solution  of  (53)  and  (54)  with  El  =  5.05. 

£(0)«  -  5.0069  +  3.26742;.  n<0)  =  -  8.631 14. n<0)  -  1.86353. and oth- 
er  parameters  as  given  in  Fig.  6. 
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FIG.  10.  Solution  of  (53)  «nd  (54)  with  El„  -  10.0,  £(0)  -  -  6.0  +  2.W. 
KO)  m  —  10,  >i<0)  »  0,  and  other  penmeten  u  given  in  Fig.  6. 


where  Eg  is  a  real  nonnegative  parameter. 

To  illustrate  the  qualitative  features  of  the  solutions  of 
the  above  equations,  numerical  results  are  obtained  for  the 
case  where  m  »  1,  L  *»  v/VlQ,  /',  «•  2.0,  and  y,  =*  1.0. 
These  values  are  chosen  to  simplify  the  numerical  computa¬ 
tion.  They  may  not  correspond  to  any  particular  physical 
situation.  Numerical  results  for  specific  physical  situations 
will  be  presented  elsewhere. 

First,  we  compute  the  locus  of  the  equilibrium  electric 
field  E’  as  a  fhnction  of  E  J  using  (27)  and  (20).  Figure  1 
shows  the  locus  in  the  (£*  *£,)-plane.  It  can  be  seen  that 
bifurcation  occurs  at  E J  «■  E&  »  0.09974.  At  E\  »  E ^ 

=»  5.05,  one  of  the  equilibrium  points  E '  returns  to  the  ori¬ 
gin.  AsEo  —  ao,  the  locus  of  one  of  the  equilibrium  points  is 
asymptotic  to  the  line  E,  =■  (2ct, 

*  -0.1  Er,  while  the  other  one  tends  to  (ErJi,) 


=  (  —  oo  ,0).  Figure  2  shows  the  equilibrium  ion  densities  n‘ 
as  a  function  of  E  l . 

Next,  we  examine  the  nature  of  each  equilibrium  state 
z,  for  various  values  of  E  a  by  determining  the  roots  of  (37) 
or  the  eigenvalues  of  Jr  (z,  ;£n ).  Figure  3  shows  the  locus  of 
the  eigenvalue,  of  J  r  (0 \Eg )  as  El  varies.  As  established  in 
Sec.  4,  when  £o  <EL  —  5.05,  all  the  eigenvalues  have  nega¬ 
tive  real  parts  implying  that  the  origin  is  asymptotically  sta¬ 
ble.  When  E  l  exceeds  E&,  one  of  the  real  eigenvalues 
crosses  the  imaginary  axis.  Consequently,  the  origin  be¬ 
comes  unstable.  Figure  4  shows  the  eigenvalue  locus  of  i  f 
(z *  \E0)  for  0. 1  <£  o  <40.0,  where  the  components  of  z* 

=  (n*  ,0,£  V Eg 4*  E0)  are  given  by  ( 1 7),  (20),  and  (27). 
We  note  that  for  0. 1  <E £  <0. 108,  all  the  eignevalues  have 
negative  real  parts,  and  at  is  J  =0.  108,  a  complex  conjugate 
pair  of  eigenvalues  cross  the  imaginary  axis  into  the  right- 
half  plane.  It  can  be  verified  that  Hopf  bifurcation  takes 
place  at  this  point.  The  locus  of  the  eigenvalues  of  Jf  (z~ ; 
Eg )  is  shown  in  Fig.  5.  Here,  for  0. 1  <.£  J  <2.669,  J  f  (z,*  JE0) 
has  a  positive  real  eigenvalue.  When  E  l  >2.67,  all  the  eigen¬ 
values  of  (z~  \Eg  )  are  in  the  left-half  plane. 

An  inspection  of  the  eigenvalue  loci  given  by  Figs.  3-5 
suggests  that  one  might  search  for  the  existence  of  chaotic 
solutions  in  the  neighborhood  of  z*  for  E\>  0. 108  (Hopf 
bifurcation  point).  Numerical  integration  of  (53)  and  (54) 
with  various  initial  conditions  was  performed  for  progres¬ 
sively  larger  values  of  E\.  The  results  suggest  that  the  peri¬ 
odic  solutions  in  the  neighborhood  of  z,+  (whose  existence  is 
ensured  by  the  Hopf  bifurcation  theorem)  are  unstable  and 
the  bifurcation  is  subcriticaL  Figure  6  shows  the  time-do- 
main  buildup  of  a  nearly  periodic  solution  which  evolves  into 
chaotic  oscillations.  The  projection  of  the  trajectory  onto  the 
(£*,2f,)-plane  is  shown  in  Fig.  7.  Figures  8-10  show  the 
chaotic  solutions  for  various  values  of  E\ .  It  was  found  that 
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(55) 


+  ^r  +  (ji'~  +  an)2  - 
=  —  n(p2  +  an)E0, 

d%  •  (2r_ 

V  (//*  4-  an)  /  i 

4  [f  +  (p2  4-  an)2  - 


aryn 


dr 


=  -(»  +  (>■- 


an 


(p2  4-  an) 


ayt  \ 
(p2  4-  an)  / 

■)njE0- 


h 


E, 


(/i2  4-  an)  j 

Considering  n  as  a  slowly  time-varying  parameter,  the  fre¬ 
quency  of  electric-field  oscillations  is  roughly  equal  to 
u  =  [yfy2  4-  an)2  —  ayh(p2  4-  an) -l  ] 1/2,  and  the  effec¬ 
tive  damping  coefficient  is  2 y  —  ah(p2  4-  an)  -l .  Let  T  de¬ 
note  the  time  interval  corresponding  to  an  ion  density  dip, 
and  t  *  is  the  minimum  point  of  n  over  T  where  ri(r  •)  =  0. 
Then  co(t  •)  >  oAt )  for  all  r  in  T  such  that  /i(r )  and 
p 2  4-  an(t )  have  the  same  sign.  This  condition  is  satisfied  for 
the  solutions  shown  here. 

Figure  12  shows  the  power  spectra  of  the  electric  field 
computed  by  means  of  the  fast  Fourier  transform  method. 
The  results  resemble  those  corresponding  to  turbulence. 
Also,  the  spectral  bandwidth  increases  with  £  £  as  expected 
from  physical  considerations.  Finally,  the  truncated  discrete 
version  of  the  autocovariance  function  of  E  given  by 


4>.00  10. 00  JO. 00  30.00  10.00  SO. 00  60.00 

FREQ.  -10' 

FIG.  1 2.  Power  spectra  of  the  electric  field  corresponding  to  the  solutions  of 
(53)  and  (54)  as  shown  in  Figs.  8-10;  frequency  scale:  1/40.96  normalized 
unit. 


PU*  )  = 


1 


£'  [E(iA)  —  E] 


(N-j+\) 

X{E[(i  +  j-\)A  ]-£}• 


(57) 


these  solutions  are  highly  sensitive  to  initial  conditions. 

Also,  not  all  trajectories  in  the  z  space  tend  to  the  chaotic 
solutions  as  r  -*  ®> .  This  is  apparent  from  the  fact  that  for 
0.108<£|  <£^  —  5.05,  the  origin  z  »  0  is  a  stable  equilib¬ 
rium  point,  and  for  any  El  >  Elc  ,z,  is  always  a  stable  equi-  ^  CONCLUDING  REMARKS 


is  computed,  where  £  denotes  the  mean-value  of  E,  and  A  is 
the  time-step  size.  Figures  13a- 13c  show  the  real  and  imagi¬ 
nary  parts  of p{jA  )  for  E\  —  1.669, 5.05  and  10.0.  It  can  be 
seen  that  both  Re  p{jA  )  and  Im  p(jA  )  decay  from  their 
maximum  values  and  then  fluctuate  about  zero.  But  we  can¬ 
not  deduce  that  the  autocovariance  function  actually  tends 
to  zero  as  the  time  delay  r  — ►  oo  as  in  the  case  of  solutions  on 
a  strange  attractor. 


librium  point.  We  note  from  Figures  8b-10b  that  in  each 
case,  there  exists  a  circle  with  minimum  radius  rmin  which 
encloses  the  projection  of  the  trajectories  onto  the  (£*  ,£,)- 
plane.  Figure  1 1  shows  the  variation  of  as  a  function 
|£0 1  as  obtained  from  the  numerical  solutions.  Evidently, 
can  be  bounded  by  a  linear  fundtion  of  |£0  | .  This  is 
consistent  with  the  estimate  $  given  in  Sec.  5.  It  can  be  read¬ 
ily  verified  that  in  each  case,  the  projection  of  the  trajectories 
of  the  chaotic  oscillations  onto  the  (£*  ,£,)-plane  is  com¬ 
pletely  contained  in  Si  ■»  {£:  |£  4-  a  £0  | 

<t/2o  *  8.3445 1£0|>,  where  a  *  1.204367.  Also,  we  ob¬ 
serve  from  these  solutions  that  the  maximum  depth  of  the 
ion  density  troughs  increases  with  £  l ,  and  the  electric  field 
oscillates  more  rapidly  during  the  ion  density  dips.  This  can 
be  roughly  explained  by  considering  the  following  equations 


for  Eh  and  Ei 

derived  from  (31): 

r‘4.  (,r 

ah  > 

i  dE„ 

~*\2r 

(p2  4- an)/ 

1  dt 

It  was  found  that  the  single-mode  equations  derived 
from  the  Zakharov’s  model  for  Langmuir  turbulence  in  a 
plasma  with  phenomenological  damping  exhibit  chaotic  so¬ 
lutions  whose  power  spectra  have  turbulence-like  features. 
In  the  case  of  multiple  modes,  if  all  the  mode  coupling  terms 
are  omitted,  then  we  obtain  sets  of  uncoupled  equations  of 
the  form  (9)  and  (10).  Each  set  is  capable  of  producing  chaot¬ 
ic  solutions  when  £  l  exceeds  a  certain  threshold  value  (gen¬ 
erally  different  for  each  mode).  The  total  power  spectrum  of 
the  electric  field  is  simply  the  sum  of  the  single-mode  power 
spectra.  This  seems  to  imply  that  energy  transfer  between 
various  modes  is  not  necessary  in  producing  turbulence 
which  is  contrary  to  the  cascade  theory  of  turbulence.  There 
are  a  number  of  computer  studies’  of  Langmuir  turbulence 
induced  by  interacting  collapsing  solitary  waves  based  on 
Zakharov’s  model  with  phenomenological  damping.  Per¬ 
haps  these  computer  results  actually  correspond  to  some 
form  of  chaotic  solutions  which  are  inherent  in  the  model. 
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In  this  work,  we  have  sought  solutions  in  terms  of  the 
eigenfunctions  of  the  Laplacian  operator  over  a  bounded 
spatial  domain.  Of  course,  we  may  expand  the  solutions  in 
terms  of  any  suitable  countable  basis  for  L  2(/2 )  and  arrive  at 
a  countably  infinite  system  of  ordinary  differential  equations 
similar  to  that  given  in  (5)  and  (6).  One  may  also  consider 
directly  the  Hopf  bifurcation  problem  for  Zakharov's  model 
(1)  and  (2)  without  resorting  to  modal  expansions.  Some 
results  in  this  direction  have  been  obtained  recently.  They 
will  be  reported  elsewhere. 

Finally,  we  note  that  the  presence  of  the  phenomeno¬ 
logical  damping  coefficients  ym  and  r„  in  the  simplified 
equations  for  each  mode  m  is  essential  for  the  existence  of 
chaotic  solutions.  But  there  does  not  exist  a  clearcut  way  of 
introducing  the  damping  terms  into  the  Zakharov’s  model 
based  on  physical  considerations.  Also,  a  detailed  study  of 
the  structure  of  the  stable  and  unstable  manifolds  associated 
with  the  equilibrium  states  is  necessary  for  revealing  the  na¬ 
ture  of  the  chaotic  oscillations  described  here.  Unfortunate¬ 
ly,  this  task  is  complicated  by  the  system's  dimensionality. 
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APPENDIX 

From  (33),  it  is  evident  that 

det[J/(2*^J)  -  A  IKA,+1)  ]  »  A  (2T  +  A  )detP  -  detQ, 

(Al) 

where 

[  -  +  an')I„  — (r-M)I,J’ 

2(0Ei -fi2E0)r  2/?(EJ)r' 


-M _ 

aEJ 

l-(E0+aEl) 


(A2) 


SincedetP=  [(y  +  /t)2  +  (M:  +art€):]‘v>0,  P~‘  ex¬ 
ists  and  is  given  by 

P  =  t(y  +  k  )2  +  (At2  +  anO2] 

—  (y  +  /i  )I  v  -  (M2  +  anOI.v  j 
X  .  (m2  +  anOI.v  -  (r  +  A  )I.v  J 
Now,  det  Q  can  be  computed  by  considering  the  matrix 

,  1  OjJ 

S  =  Q 


(A3) 


o;,  p-1 
•} 

-r  _ 
aEf, 


-(E0+aVK) 


2[(^e;  zSEo)TJ<&)r]*_ 

Izv 


(A4) 


Since  det  S  =  (det  Q)(det  P  "' ),  det  P  1  =  [{y  +  A.  )2 
+  0t2  +  an,)J]_A'and 

aE?, 


detS  =  -m:-2  [(0E'k  -M2£0)r]A-' 


L-(£0+aEi) 

(A5) 


we  have  det  Q  =  det  S/det  P  1 .  The  expression  (34)  is  ob¬ 
tained  directly  from  (Al)  and  (A5). 
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ABSTRACT 


Further  studies  are  made  on  the  chaotic  oscillations  of  the  single 
mode  equation  derived  from  Zakharov's  model  for  Langmuir  turbulence  in  a 
plasma  in  the  presence  of  ar.  external  spatially  homogeneous  electric 
field  oscillating  at  the  electron  plasma  frequency.  First,  the  bifur¬ 
cation  of  equilibrium  states  with  respect  to  both  the  energy  of  the  exter¬ 
nal  electric  field  and  the  damping  coefficient  of  the  high-frequency 
waves  is  investigated.  Then,  studies  are  made  on  the  onset  and  quench¬ 
ing  of  chaotic  oscillations  when  damping  is  varied.  The  nature  of  the 
chaotic  oscillations  is  explored  by  determining  various  Poincare  maps 
numerically. 
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I.  INTRODUCTION 


Recently,  it  was  shown^  that  the  single  mode  equation  derived  from 
2 

Zakharov's  model  for  Langmuir  turbulence  in  a  plasma  in  the  presence  of 
an  external  spatially  homogeneous  electric  field  oscillating  at  the 
electron  plasma  frequency  and  with  phenomenological  damping  has  nonperiodic 
chaotic  solutions.  The  power  spectra  corresponding  to  these  solutions 
have  turbulence-like  features.  In  reference  1,  the  chaotic  solutions 
were  discovered  through  bifurcation  analysis  and  numerical  experimentation. 
The  structure  and  the  onset  of  these  solutions  were  not  explored  in  detail. 
Here,  further  studies  are  made  on  the  nature  of  the  chaotic  oscillations. 
First,  the  bifurcation  of  equilibrium  states  with  respect  to  both  the 
energy  of  the  external  electric  field  and  the  damping  coefficient  of  the 
high-frequency  waves  is  studied.  Then,  the  quenching  of  chaotic  oscil¬ 
lations  due  to  increased  damping  is  explored.  Finally,  some  properties 
of  the  chaotic  oscillations  are  studied  by  determining  various  Poincare 
mappings  numerically. 


II.  SIMPLIFIED  MODEL 

As  developed  in  reference  1,  the  simplified  equations  for  a  single 
mode  m  in  the  Zakharov's  model  have  the  form: 


dE 
_ n 

dt 


(u  -  iy  ) E  =  n  (E  +  a  E  ), 
m  mm  mo  mm 


d2n 


dn 


dt 


+  2r  -7-^  +  y  n  =  -y  E  •  (E  +  E*)  +  3  E  ; 
2  mdt  mm  mo  m  m  m  m1 


(1) 

(2) 


with 


a 

m 


<in. 


L 


v2t$f,(xm 


m 


(x)  dn, 


(3) 


where  a-b  denotes  the  usual  scalar  product  of  two  vectors  a  and  b  in  the 


1 


N 

real  n-dimensional  Euclidean  space  ]R  ;  is  the  orthonormalized  eigen- 

m 

function  of  the  Laplacian  operator  corresponding  to  the  eigenvalue  -y^  <  0 
Y  and  f  are  the  phenomenological  damping  coefficients:  E  and  n  are 

mm  to  mm 

respectively  the  coefficients  of  expansions  for  the  electric  field  E  and 
ion  density  n: 

E(t,x)  *  ^^Ev(t)<j)lf(x),  n(t,x)  =■  y  ^  nk(t)<J)k(x)  (4) 

k  k 

N 

defined  on  the  spatial  domain  ft  C  M  .  E  *  (E  ,,...,E  „)  is  a  real 

o  ol  oN 

constant  vector  representing  the  normalized  amplitude  of  the  external 
electric  field.  In  what  follows,  we  shall  omit  the  subscript  m  in  (1) 
and  (2)  for  brevity. 

Let  Er  =*  Re(E),  E^.  -  Im(E) ,  and  z  denote  the  2 (N+l) -dimensional  real 
T  T 

vector  (njhjE^Ej)  ,  where  (•)  denotes  transposition.  Equations  (1) 
and  (2)  can  be  rewritten  as: 


dz 

dt 


f (z;y»Eq)  - 


Vn  -  2rn  -  2u2Eq.Er  +  0(  ||ER||2  +  ||EI||2) 

-yE^  +  (y2  +  an)Ex 

-YEj  -  (y2  +  an)ER  -  qEq 


where  ||  •  j|  denotes  the  Euclidean  norm. 


(5) 


III.  BIFURCATION  OF  EQUILIBRIUM  STATES 

The  equation  for  the  equilibrium  states  is  given  by 


f(z;y,Eo)  -  0,  (6) 

2  (U+l) 

where  y  and  Eq  are  parameters.  Evidently,  z«0  £  S  is  an  equili¬ 

brium  state  for  all  Y>EQ.  From  the  implicit  function  theorem,  a  nece¬ 
ssary  conditions  for  bifurcation  from  the  main  branch  C*  {(0,y,E  ):y€®» 

o 

Eq  £  R  }  is  that  the  Jacobian  matrix  Jg(0;Y,EQ)  given  by 


2 


Jf(0;y,Eo) 


0 

1 

oT 

T 

°N 

-y2 

-2r 

-2y2ET 

o 

°N 

0 

0 

-^N 

y2I 

M  N 

-E 

o 

0 

-lJlN 

-yIn 

(7) 


is  singular,  where  0N  and  I  are  the  N-dimensional  zero  vector  and  N*N 
identity  matrix  respectively. 

Since 


det  Jf(0;y,Eo)  =  (y2  +  y4)1 


“{y2(y2+ 


y4)  -  2y4||  E 


!h 


(8) 


bifurcation  from  C  can  occur  only  when  det  (0 ;y , Eo)=0  or 


||2=  E2c  £  (y2+  y4)/ (2y2) . 


(9) 


zero 


It  can  be  deduced  from  the  following  exact  expressions  for  the  non- 

e  e 

equilibrium  states  zg  =  (n^O.E^.E^)  with 

ER  =  StEo’  EI  =  ^iE0>  (10> 

g*  =  5x(y2  +  8y-1||Eo|P^)[y  +  (2a  +  By~2)||  Eq  H2?*]-1 ,  (11) 

f±  =  _ ^ _ J||E  II2  -  (m28  +  y26)  . 

1  (B2  +  62y2)||Eo!|2t  °  6y2 


(||Eo||4-  4aB6"2||Eo||2  -  (^f1)2]7},  (12) 

n*  =  (By-2 [<?*) 2  +  (5*)2]  -  2g±  }||Eo||2,  (13) 

||E0I|2  >  E2c  &  2a6-2 [3  +  (B2+62y2)2]  (14) 

6  -  (2a  +  By-2),  (15) 


that  bifurcation  from  t  indeed  takes  place  at  any  point  (0,y,EQ)£  such 

that  || E  II  *  E  .  The  critical  external  field  energy  E2  increases  with 
"  o"  oc  oc 

the  damping  coefficient  y. 

For  any  fixed  y  and  0  ^  i!E0l[2  <  E2c>  the  origin  z  =  0  is  the  only 
equilibrium  state.  When  ||Eo||2  =  E2^,  a  nonzero  equilibrium  state  zg 


3 


emerges,  and  as  II E  II2  increases  from  E2  ,  z  bifurcates  into  two  distinct 

O  OC  Q. 

equilibrium  states  z^  and  zg.  Moreover,  the  bifurcated  branch  correspond¬ 
ing  to  z*  crosses  over  the  main  branch  C  at  if2  .  Figure  1  shows  the  bi¬ 
furcation  diagram  for  E®,  the  imaginary  part  of  the  equilibrium  electric 
field,  with  j | Eq 1 12  as  a  variable  parameter  for  the  one  dimensional  case  (N»l) 
with  a  bounded  spatial  domain  Q  -  [0,L]  with  (x)  ■  (2/L)2sin  (trx/L)  and 
parameters : 

y..  *  it/L,  L  »  7T//10, 

1  1 
ax  *  (27/L)2/(3tt), 

Evidently,  the  new  branch  appears  abruptly  when  1 1  Eq  J  j2  *  E2c>  and  immediately 
bifurcates  into  two  distinct  branches  as  II E  ||2  increases  from  E2  . 

II  Q  II  QQ 

Now,  we  consider  the  case  where  1 1  Eq 1 12  is  fixed  and  the  damping  coef¬ 
ficient  y  is  a  variable  parameter.  It  can  be  readily  verified  that  the 
trivial  solution  z  =  0  is  the  only  equilibrium  point  when  y*  0.  From  (12), 
it  is  evident  that  for  real  £*,  y  must  satisfy: 

A  ||  E  ||2  i 

Y  *  Y  " - - —  (!Ie0II2  (2cty2+  S)2-  4aBy4}2  .  (17) 

C  2ay2  ° 

For  0  <  y  <  Yc»  there  exist  two  distinct  nonzero  equilibrium  states 

and  they  emerge  into  one  when  y  “  Y  •  When  y  increases  beyond  y  ,  the  non- 

C  ~  C  - 
zero  equilibrium  state  abruptly  disappears.  When  y  ■  yc  *  y £ 2( j  Eq J f2  -  y2}2, 

there  is  only  one  nonzero  equilibrium  state.  Figures  2  and  3  show  the 

0  6 

loci  of  equilibrium  E^  and  n  as  y  increases  from  zero  for  the  one  dimen¬ 
sional  case  with  parameters  given  in  (16)  and  fixed  values  of  ]|  Eq | p .  We 
observe  from  Fig. 3  that  as  y-+0,  ne  tends  to  nonzero  values  which  do  not  cor¬ 
respond  to  equilibrium  ion  densities.  Figure  4  shows  the  loci  of  equili¬ 

brium  electric  field  as  y  varies  while  keeping  Eq ||2  constant. 


Y1  "  lf  ^  (16) 
\  -  -8tt(2/Ls)z/3. 


4 


respect  to  parameter 


2.669,  5.05 
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IV.  QUENCHING  OF  CHAOTIC  OSCILLATIONS 


In  Che  Zakharov's  model  for  Langmuir  turbulence,  there  does  not  exist 

a  clearcut  way  of  introducing  damping  terms  into  the  equations  based  on 

physical  considerations.  Here,  phenomenenological  damping  coefficient  ym 

and  are  introduced  into  the  simplified  model  (l)-(3).  They  may  corre- 

3 

spond  to  Landau  damping  of  the  high  and  low  frequency  waves  given  by 


£)(f)\ 


A.) 


-3 


exP[-  j(MmADr2  -  f], 


r  _  2  /  h\t 
m  ~  4 


(18) 


(I)Vd* 


where  is  the  Debye  length  and  m^/m^  is  the  ion-electron  mass  ratio.  We 

observe  that  the  damping  coefficient  Y  varies  rapidly  for  small  y  L  >  0. 

m  m  u 

Therefore,  it  is  of  interest  to  determine  any  change  of  the  behavior  of 

the  solutions  as  Ym  is  varied,  in  particular,  the  onset  or  quenching  of 

chaotic  oscillations  when  y  tends  to  some  threshold  value.  This  is  studied 

'm 

numerically  for  the  one-dimensional  case  with  the  fixed  parameters  as  given 
in  (16). 

First,  we  study  the  case  with  E2  =  2.669.  It  is  known  that  chaotic 
oscillations  exist  for  V  ■  2.0  and  y=  1.0  (see  Ref.l).  At  these  para¬ 
meter  values,  the  equilibrium  states  consist  of  the  origin  z  =  0  and  two 
distinct  nonzero  equilibrium  states  z^  and  zg  given  by  (10)-(15).  Chaotic 
oscillations  were  found  in  some  neighborhoood  of  zg.  Figures  5-7  show 
the  loci  of  eigenvalues  of  J^z^.YjE^  with  »  2.669  and  variable  para¬ 
meter  y  for  each  of  the  equilibrium  states  zg.  We  observe  that  at  y ■  1.0, 
there  exist  a  pair  of  unstable  complex  eigenvalues  for  J f(ze»Y»E0)  with 
E2  ■  2.669.  As  Y  is  increased,  these  eigenvalues  cross  the  imaginary  axis 
at  ymA.35.  When  y  is  increased  beyond  4.35,  all  the  eigenvalues  of 

4*  <% 

Jf (z^jYjEq) ,  E*  =  2.669,  remain  in  the  left-half  complex  plane.  By  defin- 
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Fig. 5  Loci  of  eigenvalues  of  Jf(0,y,Eo)  with  IT  =  2.669,  and  variable 
parameter  y. 
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ing  a  parameter  y  =  1/y,  we  see  that  Hopf  bifurcation  with  respect  to 

parameter  y  takes  place  at  y  ~  1/4.35.  When  y=  y  =  5.3358,  one  of  the 

c 

real  eigenvaues  is  zero.  At  this  point,  z^  and  z^  coalesce  into  one  equi¬ 
librium  state  as  shown  in  Figs.  2  and  3.  Thus,  we  expect  that  quenching 

of  the  chaotic  oscillations  will  occur  when  y*4.35  for  E2  =  2.669.  Nu- 

o 

merical  integration  of  the  system  equations  confirms  our  expectation. 

Figures  8a  and  8b  show  the  existence  of  a  periodic  solution  at  y=  4.35  as 
predicted  by  the  Hopf  bifurcation  theorem.  As  the  damping  y  is  reduced, 
the  amplitude  of  the  periodic  oscillations  increases  as  shown  in  Figs.  9a 
and  9b  with  y*  4.30.  At  y=  3.0,  a  more  complex  form  of  periodic  solutions 
appears  (Figs. 10a  and  10b).  As  y  is  further  reduced,  the  oscillations 
become  more  complex  in  structure  and  appear  to  be  almost  periodic  for 
2.0  >  y  >  1.6  (see  Figs. 11-13).  Finally,  for  y<1.55,  the  solutions  be¬ 
come  chaotic  (see  Figs.  14a  and  14b).  Figures  15-17  show  the  chaotic  oscil¬ 
lations  for  y=  1.5, 1.4  and  1.1.  The  power  spectra  of  the  electric  field 
corresponding  to  the  solutions  in  Figures  14-17  are  shown  in  Fig. 18.  It 
can  be  seen  that  the  spectra  evolve  from  essentially  discrete  spectra  to 
broadband  spectra  as  y  crosses  the  threshold  value  YT  ~  1-35.  Although 
the  solutions  shown  in  Figs. 10-17  start  with  the  same  initial  conditions, 
similar  behavior  has  been  observed  for  solutions  starting  near  the  equili- 

-  brium  state  z+. 

e 

V.  POINCARE  MAPS 

To  obtain  some  idea  on  the  nature  of  the  chaotic  oscillations,  we  con¬ 
sider  various  Poincare  maps  associated  with  (4)  with  parameters  given  in 
(16).  Let  H  denote  a  given  three  dimensional  hyperplane  in  the  state 

It, 

space  B.  of  system  (4)  whose  flow  is  denoted  by  the  family  of  mappings 

4  4  / 

F^  :  ]R  -*•  31  ,  t€lR.  We  define  a  Poincare  or  first  return  map  P  from 


13 


Solution  of  (5)  with  initial  condition:  E(0)  a  -3.07656  +  4.223281 
n(0)  *  -8.38056,  n(0)  -  -8.25344;  E2  -  2.669;  y  -  4.35. 


J/> 


<JT 


S  Q  H  into  S  such  that 


P(z)  =  F  .  .  (z)€S  for  all  z  G  S,  (19) 

T  \  Z  ) 

and  S  is  transverse  to  the  flow  (i.e.  (3F  (z)/3t)|  .  >  has  a  nonzero  com- 

t  'T(z) 

ponent  with  respect  to  n,  the  normal  vector  for  H) ,  where  T(z)  is  the  small¬ 
est  time  t  >  0  such  that 

Ft(z)(2>^S  and  s8nt (z-h) (F,^  v (z) • n) 1  >0.  (20) 

Let  z(0)  be  a  point  in  S  at  t  =  0,  we  are  interested  in  the  sequence 
k+1 

of  states  z  generated  by 
k+1  k 

ZK  =  P(zK),  k*0,l,2, . . . , 

0  (21) 
z  =  z (0)  G  S. 

If  there  exist  a  point  z°  and  an  integer  K  >  0  such  that  zK  =  z°,  then  z° 

K  k 

corresponds  to  a  periodic  solution  with  period  T  *  t(z  ).  For  non- 

k=0 

periodic  or  chaotic  oscillations,  there  do  not  exist  such  points. 

First,  we  obtained  numerically  the  points  generated  by  Poincare  maps 
corresponding  to  various  hyperplanes  in  the  state  space  of  system  (4)  with 
parameters  given  in  (16)  and  F  =2.0  for  various  values  of  E‘.  Figures 

19-22  show  the  points  generated  by  Poincare  maps  for  progressively  higher 
values  of  E*  and  fixed  damping  coefficients  y-  1.0  and  T=  2.0  along  chaotic 
solutions.  A  notable  feature  of  these  results  is  that  the  points  appear 

^  j 

to  lie  along  certain  curves  in  K  ,  which  implies  that  the  Poincare  maps 
associated  with  chaotic  solutions  are  one  dimensional  in  nature.  But  the 
graphs  of  the  Poincare  maps  in  terms  of  some  curve  parameter  are  not  readily 
obtainable. 

VI.  CONCLUDING  REMARKS 

The  results  of  this  study  show  more  clearly  the  onset  or  quenching  of 
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Fig. 20  Isometric  representation  of  points  generated  by  Poincare  maps 
onto  the  hyperplanes  (a)  *  0;  (b)  E^  *  0;  (c)  dn/dt  =  0. 


Angle  of  third  axis  is  at  45°;  E*  ■  2.669. 
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chaotic  oscillations  for  fixed  E2  and  variable  damping  coefficient  y  as 

2  1 

compared  to  the  case  with  fixed  y  and  variable  Eq  studied  earlier.  As 

/■* 

y  is  decreased  from  the  Hopf  bifurcation  point  y,  the  oscillatory  solutions 
evolve  from  simple  periodic  to  almost  periodic  solutions  and  to  chaotic 
solutions.  It  is  conjectured  that  the  chaotic  solutions  appear  after 

A 

further  bifurcation  from  the  periodic  solutions  as  y  decreases  from  y.  The 
verification  of  this  conjecture  requires  further  study. 

In  this  study,  only  the  single  mode  model  is  considered.  The  ex¬ 
istence  of  chaotic  oscillations  in  a  multiple-mode  model  is  being  investi¬ 
gated.  The  results  will  be  reported  in  the  near  future. 
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ABSTRACT 


A  numerical  study  is  made  on  the  chaotic  solution  of  the 
nonlinear  interaction  of  three  positive -energy  electrostatic  waves 
in  a  plasma,  where  one  wave  is  growing  linearly  and  the  remaining 
two  are  linearly  damped.  The  reduction  from  the  six-dimensional 
system  to  a  three-dimensional  system  and  bounds  for  the  chaotic 
attractor  are  also  discussed. 


I.  INTRODUCTION 


In  the  conventional  theory1’2  of  turbulence  in  plasmas,  turbulence 
is  described  as  the  state  in  which  a  large  number  of  collective  degrees  of 
freedom  are  strongly  excited.  The  energy  of  the  unstable  modes  is  distrib¬ 
uted  to  other  modes  by  a  nonlinear  process  and  then  dissipated  through  some 
form  of  damping  mechanism.  When  the  energy  transfer  between  the  modes  is 
balanced,  a  stationary  broad  spectrum  corresponding  to  stationary  turbulence 
appears.  If  many  collective  degrees  of  freedom  are  excited  simultaneous ly 
from  the  level  of  thermal  fluctuation,  the  random  character  of  the  fluctua¬ 
tion  is  preserved  to  some  extent.  If  only  a  small  number  of  degrees  of 
freedom  are  excited  in  the  initial  stage,  there  must  be  a  randomization  pro¬ 
cess  by  which  the  system  goes  into  the  turbulent  state. 

Recently,  it  has  been  proposed3'5  that  the  turbulent  state  or  chaotic 
behavior  can  arise  in  deterministic  nonlinear  models  which  have  a  small 
number  of  degrees  of  freedom.  Particular  examples  of  such  models  have  been 
found  in  plasma  physics6'8.  In  this  paper,  we  consider  a  simplified  model 
governing  the  nonlinear  resonant  interaction  of  three  waves  with  both  linear 
growth  and  damping  rates: 


dAj/dt  =  y1Al  -  iVAjjA, 
dA2/dt  =  -  Y2A2  -  iVA1Aj*  • 
dAj/dt  =  -  Y3A3  -  iVAjA.,* 


(1.1) 


where  Aj's  are  normalized  complex  amplitudes  of  waves,  V  is  real  ant! 
Yj's  are  positive  numbers.  System  (1.1),  without  the  growth  or  damping 
terms,  has  been  thoroughly  investigated.  We  shall  show  that,  under  certain 


2 


conditions  on  Yj  »  the  linear  instability  of  the  first  wave  could  lead  to 
a  turbulence- like  state  without  any  additional  assumption  for  randomization. 
A  particular  example  in  the  form  of  system  (1.1)  will  be  studied  in  detail 
numerically. 


II.  DERIVATION  OF  BASIC  EQUATIONS 

We  first  consider  three  nonlinearly  interacting  oscillators  whose 
equations  of  motion  are  given  by 


a3  -  2Y3aj  ^  uq  a!  ^  va2a3 


a2  +  2y2a2  +  up  2a2  =  va3a3  f  , 


a3  +  2y3a3  +  u>32a3  *  va3a2 


(2.1) 


where  v,  Yj’s  and  Wj's  are  real  positive  numbers.  We  assume  that  the 
linear  damping  on  growth  rates  is  small  (i.e.,  Yj  <<:uw)  and  the  nonlinear 
coupling  is  small  (i.e.,  |a3  jv  <<  w22  ,ui32)  .  Then  we  can  express  a^  in 
the  form 


aj  (t)  =  (Aj  (t)exp(iujjt)  +  A^*  (t)exp  (-icjj  t) }  (v/u. )  ^  ,  (2.2) 

where  A^  is  a  slowly  time-varying  component.  If  the  frequencies  satisfy 
the  resonance  condition  up  =  to2  +  u)3  ,  then  substituting  (2.2)  into  (2.1) 
leads  to  equations  of  the  form  (1.1),  with  V  =  (v3/4uqu)2u>3)  2  . 

We  can  also  derive  equations  in  the  form  of  (1.1)  for  the  nonlinear 
wave-wave  interaction  in  a  one-dimensional  plasma  involving  two  plasma  waves 


(denoted  by  subscripts  1  and  3"!  and  an  ion-acoustic  wave  (denoted  by 
subscript  2).  We  assume  that  the  resonance  conditions  for  hoth  the  wave 
numbers  (i.e.,  kt  = kj  +  k3)  and  the  frequencies  are  satisfied.  Here,  the 


pair  (ojj.kj)  satisfies  the  linear  dispersion  relation  associated  with  the 
j-th  wave.  We  assume  that  the  first  plasma  wave  is  excited  by  some  external 
force  and  has  a  linear  growth  rate  Yj  .  The  second  plasma  wave  and  the 
ion-acoustic  wave  are  excited  by  the  first  plasma  wave  through  nonlinear 
coupling  and  have  linear  damping  rates  y2  and  y,  ,  respectively.  We 
also  assume  that  y^  «  wu  and  that  the  nonlinear  coupling  of  wa\es  produces 
small  shifts  in  frequency  from  the  linear  values.  Moreover,  the  resonant 
interaction  between  waves  and  particles  is  negligible.  Then,  the  equations 
for  the  interacting  waves  are  described  by  (1.1).  The  amplitudes  A.  and 
the  parameter  V  are  defined  as  follows9: 


Aj  (t)  = 


V) 

~w\ 


3 Re  g(lj(k) 
3  w 


ki 4 


.  (2) , 


V  -  -  e^Ckj.k^k,), 


*4  ft)  .  j  =  1.2,3, 

i  J 


k^^dRec^]fk).  1,1/2 


(2.3) 


!.4) 


Here,  the  parameter  V  is  real,  is  the  slowly  varying  component  of 

*  he  “lectrostatic  potentional  of  the  j-th  wave,  and  e  ^  and 
•"  first  and  second  order  dielectric  constants  obtained  from  the 
■■■•;  Vlasov- Poisson  equations. 

r in r  ro  ohtam  explicit  expressions  for  the  variables  and  para- 
•  .  ier  *he  following  simple  model:  weak  electron  beams  are 
-  .  •  ; ;  ies  of  an  infinite  plasma  slab  with  thickness  L  . 
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We  assume  that  there  exist  two  eigen-modes  of  the  plasma  waves  and  one 
eigen-mode  of  the  ion-acoustic  wave,  and  that  the  velocity  of  the  beam  in 
the  positive  (resp.,  negative)  direction  is  between  the  phase  velocities 
of  the  two  plasma  waves  propagating  in  the  positive  (resp.,  negative)  direc¬ 
tion.  We  also  assume  that  the  interaction  between  the  beam  and  the  fast 
plasma  wave  can  be  neglected.  Then,  it  is  sufficient  to  consider  the  system 
consisting  of  the  beam  in  the  positive  direction,  the  plasma  wave  '1'  and 
the  ion-acoustic  wave  '2'  traveling  in  the  positive  direction  and  the 
plasma  wave  '3'  traveling  in  the  negative  direction. 

We  further  assume  that:  (i)  the  electrons  and  ions  in  the  plasma  and 
the  electrons  in  the  beam  have  time -independent  Maxwell  velocity  distribu¬ 
tions  in  the  zero-th  order,  (ii)  the  Landau  damping  of  wave  '2'  is  negli¬ 
gible  as  compared  to  the  damping  due  to  Coulomb  collisions,  but  the  Landau 
growth  rate  of  wave  '1'  due  to  the  beam  is  larger  than  the  collision  damp¬ 
ing  rate,  and  (iii)  the  plasma  is  nonisothermal  so  that  the  Landau  damning 
of  the  ion-acoustic  wave  is  negligible.  Then,  the  amplitudes  ,  the 
parameter  V  and  the  damping  on  growth  rates  Yj  are  expressed  explicitly 
as  follows: 


A.  = 
J 


k.2 

8tt 


2 

^e 


y2 


» j =  i » 3, 


,  0)  , 

_  J_  pe  _c_  kj 


8tt  v.-  2  in  a)  c 
Tc  c  pe  s 


k 2 

2u>  ■ 

A2  - 

c 

8t r 

k^c 

pj 

.X 

lL 

m  )] 
c  , 

|4TTkzCc 

m.  j 
ij 

J 

% 


(2.5) 


Vz 


(2.6) 


-  Y  7  +  /ii  aj 
col  pe 


V  v/k0 

r>Ai> 

2 

'vd  -  “pe/k>l 

l 

\  j 

\  j 

exp  ■ 

VT 

Tb 

(2.7) 
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Y2  = 


m 

tt  _e 

F  m. 


k2Cs  *  YS 


'  col 


(2.3) 


Here,  Ycoi  is  the  collisional  damping  rate,  mg  (resp.,  nr)  is  the 

electron  (resp.,  ion)  mass,  e  is  the  electron  charge,  Wpe  (resp.,  w 

is  the  electron  (resp.,  ion,  beam)  plasma  frequency,  v_  (resp.,  vT  )  is 

je  ‘b 

the  electron  thermal  velocity  of  plasma  (resp. ,  beam)  ,  cg  is  the  sound 
velocity  of  plasma  and  v^  is  the  beam  velocity. 

It  can  be  shown  that,  if  Aj  =  A2  =  0  at  the  initial  time,  then 
|At  (t)  |  -®  as  t-®  .  Hence,  the  solution  of  (1.1)  is  not  uniformly  hounded 
with  respect  to  the  set  of  initial  conditions  whose  closure  includes  the 
set  such  that  Al  =  A2  =  0  ,  Thus,  the  assumptions  made  in  the  derivation 
of  (1.1)  may  be  violated  in  some  cases.  Such  a  situation  will  be  discussed 
in  Section  4. 


III.  ANALYSIS  OF  THE  SIMPLIFIED  MODEL 

We  eliminate  the  parameter  V  from  (1.1)  by  setting  * -VA^ (j=l,2,3) 


Xj  =  Y1X1  +  iX2X3 


X2  =  -  y2X2  +  iXjX3 


X3  =  -  y3X3  +  iXjX2 


(3.1) 


The  time  derivative  of  £  =  ^{2XiXj*  +  X2X2*  +  X3X3*}  is  given  by 
K  =  2YiXiXi*  -  y2X2X2*  -  y3X3X3*  .  Hence, 


i(t)  <,  £(0)exp(2  ma»{Yi,Y2.73>|t|]  ,  te  (-«,“) 


(3.2) 


Therefore,  the  integral  curve  of  (3.2)  exists  for  all  te 

i9  • 

Setting  Xj  =  r^e  J(j=],2,3)  ,  with  positive  r^  , 


(-°°,°°)  • 

(3.1)  becomes 


i'i  +  i©i x*x  *  YjTj  +  ir2r3e  1  i ~92 ~e 3 ^ 

r2  +  ie2r2  =  -  72r2  +  ir j^e* *-01  *02 "<53-*  •  .  (3.3) 

r3  +  i03r3  =  -  y3r3  +  ir!  r2ex  ('01  ~*2  ‘03-)  J 

Taking  the  time  derivative  o^  n  =  Re(X1X2:*X3*)  =  r^r^os (9 3 -02 -e3)  , 
we  have 


n  =  -  (r2  +  r3  -  Y!)n  .  (3.4) 

Hence,  if  Yi<Y2+Y3  ,  n(t)-*0  as  t-*-®  . 

Let  Y  =  (Xi,X2,X3)  ,  then  the  set  E  =  { Y  :  n  =  0  }  is  an  invariant  set 
in  the  sense  that,  if  Y(0)  eE  ,  then  its  corresponding  full  orbit  is  con¬ 
tained  in  Z  ;  i.e.,  Y(t)  e  E  for  all  te  (-~,®)  .  The  subset 
E0  =  (Y  :  Tj  *  r^  =  0,  j^k,  j,ke  {1,2,3}}  of  E  is  also  an  invariant  set. 

We  shall  show  that,  if  yz  f  Ys  ,  the  subset  Ey  =  (Y:  eos(9i -02-93) 

=  O,s^(0i-92-93)  *  v}-E0  ,  v  =  l,-l  ,  of  E  is  an  invariant  set,  in 
this  case,  in  the  sense  that,  if  Y(0)  e  Ev  ,  then  Y(t)  e  E^  for  almost 
all  te  (-®,«)  .  Assume  that  Y(0)  e  E  .  Since  (3.1)  has  a  unique  solution 
for  each  initial  condition,  Y(t)  e  -  E0  for  all  te  (-®,®)  .  Hence,  at  any  t, 
at  most  one  of  the  rj(t)’s  is  zero.  Furthermore,  all  of  the  r^.(t)'s  are 
non-zero  for  almost  all  te  (-®,»).  For,  if  y2  i  Y3  >  Eq.(3.1)  does  not 
have  a  solution  such  that  one  of  the  X^(t)'s  is  identically  zero  on  some 
interval  in  -E0  .  Since  n(t)*0  for  all  te  (-®,®)  ,  it  follows  that 


mr,  (9j  (t)  -  92  (t)  -  93  (t))  =  0  for  almost  all  t.  (-»,■»)  .  Moreover, 

.zin  (9t  (t)  -  02  (t)  -  9 3  (t) )  =  j  for  almost  all  te  (-'»,'*>)  .  For,  if 

sir  (0!  (t)  -  02  (t)  -  03  (t))  has  a  discontinuity  at  t,  ,  then  r.  ft  Jr,  (t2)  =  0, 

J  K 

j  +  k,  j,ke  {1,2,3}  ,  since  the  real  parts  of  the  right-hand  sides  of  Eq. 

(3.3)  are  continuous.  This  implies  that  Y(tj)  e  ZQ  ,  which  is  a  contradic¬ 
tion.  Hence,  Zy  is  an  invariant  set,  in  the  sense  defined  above. 

In  what  follows,  we  consider  the  behavior  of  Y(t)  with  Y(0)  e  I,  , 

where  v =  1  or  -1.  For  t  such  that  Y(t)  e Zy  ,  (3.3)  reduces  to 


ii  + i©! r!  =  Yiti  +  vr2r3 
r2  + i92r2  =  -  y2r2  -  vr3r3  • 
i-J  +  i63r3  =  '  Y3r3  '  ^lT2 


(3.5) 


Since  ||Xj(t)||  is  finite  at  finite  t  ,  the  integral  curve  of  (3.5) 

for  t  such  that  Y(t)  e  Z^  is  identical  to  that  of  (3.4)  for  te  (-°°,»). 

Hence,  we  can  assume  that  Eq.(3.5)  holds  for  te  (-00,00)  . 

As  mentioned  earlier,  all  of  the  rj(t)'s  are  non-zero  for  almost 

all  te  Hence,  § .  (t)  =  0  for  almost  all  t  except  when  r^  (t)  =  0  . 

Since  (t)  is  continuously  differentiable  at  all  t  ,  (t)  =  9j (0)  +  (t) 

Here,  n.  (t)  is  an  integer-valued  step  function  which  has  discontinuities 

-  i0  t  CO)  irn- (t) 

when  r^  (t)  =  0  .  Let  Xj  (t)  =  vX^e  J  =  vr^e  J  ,  j  =  1,2,3  .  Then, 

Xj  is  real  and  (3.5)  reduces  to 

Xi  =  YjXj  -  x2x3 

X2  =  -  Y2X2  +  XiXj 

x3  =  -  Y3X3  +  x2x2 


(3.6) 
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The  assumption  that  y,  "  y2  +y,  also  leads  to 


ii  +  +  ii  - 

3x  3y  32  '  ' 


(y2  +  Y3  -  Yi)  <  o  , 


H.7) 


which  means  that  the  phase  volume  shrinks  uniformly.  The  Lorenz  system 
also  has  this  property3. 

The  equilibrium  points  of  system  (3.6)  are  easily  obtainable.  If 
Yi  <  0  ,  the  origin  P0  is  the  uniaue  equilibrium  point  and  it  is  asymp¬ 
totically  stable  in  the  large.  If  Yi  >  0  >  the  equilibrium  points  are 
P0  and  P^u  =  (aCYzYa)^,  u(YiYa)  ^ ,  v(yiY2)/2),  (o.U.v)  e  15=  {(1,1,1), 
(1,-1,-1),(-1,-1,1),(-1,1,-1)>  .  The  point  P0  is  a  stable  node  in  the 
(y,z) -plane  and  a  saddle  point  in  a  plane  containing  the  x-axis.  At 
Pv  ,  the  characteristic  polynomial  of  the  linearized  vector  field  of 

0  y  U 

(3.6)  is  f(A)  =  A3+A2(y2+Y3  -yi)  +  4yiY2Y3  =  0  .  Since  the  polynomial 
has  one  negative  root  and  a  pair  of  complex  roots,  f(A)  = 

(X  +  c)  {A  -  (a  ♦  ib) } { A  -  (a-ib)}  where  c>0  .  By  comparing  the  coefficients 
of  polynomials,  we  have 

2a  -  c  =  -  (y2  +Y3  ■  Yi)  t  a2  +  b2  =  2ac  ,  (3.8a) 

c(a2+b2)  =  4y  3  y2Y  3  .  (3.8b) 

Eq.(3.8b)  implies  that  a>0  .  Hence,  at  each  ,  (o,y,v)  e  15  ,  there 

o 

is  a  bistable  plane  in  which  is  an  unstable  focus. 

a,u 

Eq. (3.6)  is  invariant  to  the  interchange  of  subscripts  2  and  3. 
Hence,  without  loss  of  generality,  we  assume  that  y2  <  Y3  •  Let  X*  (•''•.y.z) 

=  (xltx2,x3).  The  time  derivative  of  a(t)  =  y(t)/z(t)  is  given  by 
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«  *  (y 3  ~  Yz)ci  +  x(l  -  a2)  ,  (5.9) 

which  is  equal  to  u(y3-y2)  at  y  =  yz  ,  u  =  l,-l  .  Since  y2  <  Y,  ,  the 

set  F  =  {X  :  |yj  >  | z | }  is  a  positive  invariant  set;  namely,  if  X(t)  «  r 

at  some  t  ,  then  X(t)  e  r  for  all  t  >  t  . 

We  shall  show  that  F  is  an  attractor  in  the  sense  that,  for  any 

initial  point  X0  ,  X(t)  converges  to  r  as  t-*®  fi.e., 

lim  inf  ||X(t)  -  X||  =  0)  .  Since  r  is  a  positive  invariant  set,  it  is 
t-*oo  Xe  r 

enough  to  show  that,  if  X(t)  «  -r  for  all  te  [0,®)  ,  X(t)  converges  to 
F  as  t-*-®  .  Let  B(t)  =  y2(t)  -  z2(t)  ,  then  the  time  derivative  of  S 
is  given  by 

8  =  2(y322  -  Y2y2)  *  2 Cy3  -  Yz)z2  -  2y23  .  (3.10) 

If  X(t)  e -r  for  all  te  [0,®)  ,  then  |y(t)  |  <  I  z(t)  |  (or  S(t)  <  0)  and 

!y(t)|  <  (Yj/Yz)  2|z(t)|  (or  0(t)  > 0)  for  all  te  [0,®)  .  Hence,  S(t)~0 

as  t-*®  .  For,  if  there  is  a  6  <0  such  that  6(t)*6  as  t-*-®  , 

then  &(t)  > -ZyzB^  >  0  from  (3.10),  since  -B(t)sz2(t)  .  This  implies 

that  B(t)  does  not  converge  as  t-*-®  ,  which  is  a  contradiction.  Thus, 

X(t)  converges  to  3F  c  r  as  t-*®  .  Hence,  F  is  an  attractor. 

If  X(t)  g  F  for  some  t  ,  the  X- trajectory  is  trapped  in 

r+  =  F  n  {X  :  y  >  0}  or  r_  =  V  n  {X  :  y  <  0}  for  t  >  t  .  Only  this  case 

was  observed  in  the  computer  experiment  and  some  numerical  results  for 

such  a  trajectory  will  be  presented  in  the  next  section.  The  set  F+ 

i  -l 

(resp. ,  rj  has  three  equilibrium  points  Po  ,  P3 j s  ,  and  P  (resp., 

i  -i 

P0  ,  P-i,-i  and  Pi,-i)  a^d  The  combination  of  these  points  is  identical 
to  that  of  the  standard  Lorenz  attractor3. 
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We  shall  determine  the  behavior  of  the  trajectories  such  that 
X(t)  e  -  T  for  all  tc-  [0,°°).  Since  Eq.(3.6)  is  invariant  to  changing  the 
signs  of  any  two  of  the  three  variables,  it  is  sufficient  to  consider  the 
set  such  that  z  >  0  .  Let  W  =  (X  :  |ax|  =ox,  |yy|  =  yy,  z  >  0}  -  W  , 

O  )  M 

a,u«  {-1,0,1}  .  Here,  W*  V _ J  (P^  ,  }  u  (x,y,z-axis)  .  The  set  W  is 

(a,u,v)en  CT’U 

the  invariant  set  and  the  X-trajectorv  with  X0  i  W  does  not  hit  W  in 
finite  time.  In  what  follows,  the  trajectory  in  W  will  be  disregarded. 
Consider  the  following  sets  (see  Figure  1): 


S  =  {X  :  (i  =  0}  n  (W0>Q  -  r) 

Wj  =  (X  :  a>0}  n  (W0>0  -  D  • 
W2  =  (X:a<0}  n  (W0>0  -  T) 


(3.11) 


We  assume  that  X(t„)  «  SnWlf0  at  some  t„  .  In  ^  ,  (o,u)  = 

(1,-1), (-1,1)  , 

I  x  |  =  ax  ,  |y]  =  ay  ,  z<0  .  (3.12) 


Hence,  the  X- trajectory  is  transverse  to  S  and  enters 
In  Wi  nW^.j  ,  -l<a<0  ,  a>0  ,  x>0  and  x>0  . 

a  increases.  Thus,  the  X- trajectory  enters  ff,  nW^, 
Let  X(Tj)  e  (x,z)-plane  n  W1>0  and  X(Tj  +  e)  e  Wj  nWul 
X(r)eW1nW1>1  for  all  Te  [ti,!,  +e]  .  In  Wj  nW,^  , 
and  x>0  .  Hence,  from  (3.9)  if  X(t)eW1nWljl  for 


• 

Hence,  from  (3.9), 
in  finite  time. 

,  where 
0<a<l  ,  a  >  0 
all  T£  [tj  +  c,t] , 


a(t)  >  (y3-Y2)a(t)  >  (y3  '  Y2MT1  +  e)  >  0  • 


(3.13) 
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This  implies  that  the  X-trajectory  enters 
it  is  transverse  to  the  (y, z) -plane.  Let 
X(t2  +  e)  e  W,  n  W.,  fl  where  X(x)  e  Wt  nW.^j 
wi  n  W- 1 ,  j »  x<0»  x  <  0  ,  y  >  0  ,  y  <  0  , 

X(x)  c  Wj  n  W.j  for  all  Te[t2+e,t]  , 


Wj  n  W_  j  j  in  finite  time,  since 
X(t2)  e  (y,z) -plane  MV0  and 
for  all  t  f  [t2  ,t2  +  r. ] .  In 
z  >  0  and  z  <  0  .  Hence ,  i  f 


x(t)  <  Yix(t2  +  e)  <  0 

y(t)  s  -  y2y(t) 

Z(t)  ^  -  Y3 2  (t) 


(3.14) 


Consequently,  if  X(t)fiWj  nW.1(1  for  all  te  [0,°°),  then  x(t)  -*■  -°°  , 

y(t)  -»•  0  and  z(t)  -*■  0  as  t  -*■«  .  Otherwise,  the  X-trajectory  hits 

SnW.j  0  in  finite  time.  Similarly,  if  X(x0)  €  S  r»  W.j  >0  at  some  r0  , 

the  X-trajectory  converges  monotonically  to  the  x-axis  (  x(t)-+°°)  as 
t-*-°°  or  hits  S  n  W,  „  in  finite  time.  Thus,  if  X(t)eW  -f  for 

all  t<  [0,“)  ,  then  as  t-*-»  ,  the  X-trajectory  converges  to  the  x-axis 
monotonically  after  some  oscillation  about  the  z-axis,  or  oscillates  about 
the  z-axis  for  all  te  [0,®)  (i.e.,  for  any  large  t  ,  there  exist  tj 

and  t2  such  that  t2  >  tl  > t  ,  X(tj) ,X(t2)  e  (x,y) -plane  and  the  X-tra¬ 

jectory  for  t«  [t,  ,t2]  circles  the  z-axis). 

Let  us  see  roughly  how  the  trajectory  oscillating  about  the  z-axis 
for  all  t  e  [0,»)  behaves.  Integrating  (3.10), 


•os 

0  '  L  z,(t)dt  ■ 


,co 


-6(0)  +  2y2 


8(t)dt  < «  . 


(3.15) 


The  second  inequality  holds  since  8(t)  <  0  .  Likewise, 


y2 (t)dt <  ® 
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By  studying  carefully  the  vector  field  of  (5.6),  we  know  that  y  (t) 

(resp. ,  2  (t)  )  takes  a  larger  value  than  (y,y3)  ^  (resp. ,  (yiY?)^2) 
at  each  time  when  the  X- trajectory  encircles  the  z-axis.  ilencc, 
y2(t)  AO  and  z2  (t)  AO  as  t-»®  .  Therefore,  for  any  small  c  >  0, 
the  time  duration  for  which  y2(t)  >e  (resp.,  z2(t)  > e)  converges  to 
zero  as  t-*-®  .  Thus,  y2(t)  (resp.,  z2(t)  )  behaves  as  a  train  of 
pulses  whose  heights  are  larger  than  Y1Y3  (resp.,  YiY2  )  and 
whose  widths  converge  to  zero  as  t  -*■  ®  . 

The  X-trajectory,  which  is  in  W1ntV.l  v  or  W2nW,  for  all 
t  €  [tn®)  for  some  ti  and  converges  to  the  x-axis  monotonically,  is 
unstable  in  the  sense  that  any  small  perturbation  of  a  can  shift  X  into 
T  or  cause  X  to  oscillate  about  the  z-axis.  For  a  =  u(y3  -Yi), 

()y|  =uy)  on  3T  and  the  X-traje ctory  is  transverse  to  S  ,  and  3T  and 
S  approach  each  other  as  jxj  -*■«  .  The  X-trajectory  which  oscillates  about 
the  z-axis  for  all  te  [0,®)  is  also  unstable  in  the  sense  that  X  is 
shifted  into  r  by  any  small  perturbation  of  a  .  For  a  =  u(y3  -y2)  , 

( | y {  =  uy)  on  3T  and  X(t)  converges  to  9T  as  t-*-®  .  Hence,  from  the 
practical  standpoint,  the  X-trajectory  which  remains  in  -  (W  u  T )  for  all 
te  [0,®)  can  be  disregarded.  In  fact,  such  trajectoiies  were  not  observed 
in  the  computer  experiment. 
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TV.  numerical  experiment 

In  the  numerical  experiment,  we  used  the  system: 


x  =  Ylx  -  yz 

y  =  -  Y2y  +  zx 

z  =  -  z  +  xy 

which  is  obtained  from  (3.6)  by  the  substitution: 


(4.1) 


y3t-t  ,  y./y3  *  y-  ,  for  j-  1,2,3 

Xj /y 3  x  ,  x2/Yj  -  y  ,  x,/Yj  -*•  z 


(4.2) 


We  set  plasma  density  n0  *  1.62  *  107  (/cm3)  ,  beam  density  n^  =  1 . 85  *  ID1* 
(/cm3)  ,  plasma  electron  temperature  Tg  *  0.25  (eV) ,  beam  electron  temp¬ 
erature  =  0.033  (eV)  ,  mu  =  1. 16  x  10  2  3  (g) ,  (Li  ion),  L=100  (cm)  , 

27rL/kj  »  7  ,  2mL/k2  =  13  and  2ttL/  | k3 1  =  6  .  Then  y.  «  u.  in  (3.6)  and 

y2  =  0.4  in  (4.1).  For  yx  satisfying  0<y1<Y2+Y3  ,  V^=  5.57x10® 

(cm/sec)  (88.2eV). 

The  trajectories  for  various  values  of  yx  e  [0.2, 0.6]  are  shown  in 
Figure  2.  The  minimum  value  of  y2  +Y3  _Yi  is  0.8  for  yi  e  [0,0.6]. 

Hence,  from  (3.7),  the  trajectories  in  Figure  2  are  on  the  set  whose  phase 
volume  is  almost  zero  and  which  resembles  a  surface. 

In  order  to  understand  qualitatively  the  dependence  of  the  trajectory 
behavior  on  Yl  »  we  use  a  simple  model  shown  in  Figure  3.  This  model  con¬ 
sists  of  two  planes  K3  and  K2  containing  the  unstable  foci  F:  and  F2  , 
respectively.  The  trajectories  spiral  slowly  away  from  F,  (resp. ,  F2  )  on 
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Kj  (resp.,  K2  )  and  jump  from  K[  (resp.,  K,  )  to  the  line  L,  (resp. ,  L,). 
Let  <  be  the  ratio  of  the  frequency  and  the  growth  rate  of  the  spiral 
trajectory.  In  the  actual  system,  <=  jbj/a  at  Pv  t  ,  (a,y ,v)  ,  and 

O  t  u 

from  (3.8), 


(<2  -  3)  3/(k2  +  1)  =  (y 2  +  Y  3  '  Y i ) 3 /'^  1  Y 2 Y  3 


(4.3a) 


ck2  =.(k2  -  3)2(k2 +9)(Y2 +Y3-Yi)2(2yi  +Yi  +Y3)  .  n 

377  T(T^^)y12Y2Y3 -  ° 


(4.3b) 


Hence,  as  Yi  increases  from  0.2  to  0.6  ,  i c  decreases  from  5.56  to 

i  -i 

3  monotonically.  Furthermore,  P1}1  and  P.1(1  move  away  from  the  x-axis 
perpendicularly  to  the  x-axis.  Taking  into  account  these  properties  of 
the  actual  system,  we  consider  five  cases  as  shown  in  Figure  3.  Here,  as 
the  system  changes  from  (a)  to  (d)  ,  <  increases  and  Fj  and  F2 
move  upward.  Thus,  the  transition  from  (a)  to  (c)  in  Figure  2  is  simi¬ 
lar  to  that  from  (a)  to  (c)  in  Figure  3. 

In  Figure  3(a)  ,  the  focus  Ft  (resp.,  F2  )  is  below  Lj  (resp.,  l2  ). 
Therefore,  there  may  exist  various  types  of  trajectories  shown  in  the  figure 
and  their  combination  leads  to  the  trajectory  as  shown  in  Figure  2(a). 

Since  the  trajectory  can  be  arbitrarily  close  to  the  origin,  the  attractor 
may  be  unbounded.  In  Figure  3(b),  F:  (resp.,  F2  )  is  above  L!  (resp. ,  L2). 
Consider  the  trajectories  tangent  to  Lj  and  L2  at  Gj  and  G2  ,  re¬ 
spectively  and  the  shaded  region  which  they  surround.  This  region  is  ob¬ 
viously  a  local  attractor  in  the  sense  that  any  trajectory  sufficiently 
close  to  the  region  enters  the  region  at  some  time  and  remains  there. 
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This  case  corresponds  to  that  in  Figure  2(b).  As  <  takes  on  larger  values, 
the  trajectory  tangent  to  Lj  (resp. ,  L,)  becomes  closer  to  G,  fresp.,GJ. 
.Assume  that,  for  some  value  of  <  ,  the  point  B,  is  above  A2  ,  as  shown 
in  Figure  3(c).  Then,  the  local  attractor  of  type  (b)  is  replaced  by 
that  of  type  (c).  In  Figure  3(c),  we  have  omitted  another  local  attractor 
which  is  synmetric  to  the  one  shown  under  lSO°-rotation  about  the  v-axis. 

We  have  also  omitted  the  points  ,  D, ,  and  G;  on  Lx  and  E,  and 

H2  on  L2  ,  which  are  in  the  omitted  local  attractor.  If  the  points  B. 
and  C.  are  located  between  lb  and  lb ,  j  =  l,2  ,  then  these  two  attrac¬ 
tors  are  separated  and  interlinked  with  each  other.  This  case  corresponds 
to  Figure  2(c).  Assuming  that  k  is  sufficiently  large  so  that  D2  coin¬ 
cides  with  G2  as  shown  in  Figure  3(d),  the  local  attractor  of  type  (c) 
disappears. 

In  Figure  2(d)-2(f),  stable  interlinked  double  limit-cycles  and  a 
single  limit-cycle  are  observed.  They  have  a  small  numeber  of  loops.  In 
order  to  understand  their  formation,  we  consider  the  one-dimensional  re¬ 
turn  maps  of  trajectories  of  the  actual  system.  Figure  4(a)  shows  the  graph 
of  the  value  x  of  x  on  the  plane  T+  =  (X  |  z=0  ,  x,y  >  0}  with  respect 
to  the  previous  value  of  x  on  the  same  plane.  In  Figure  4(a),  the  circled 
points  indicate  the  positions  of  the  observed  limit  cycles.  Actually,  as 
shown  in  Figure  4(b),  the  mapping  from  a  sufficiently  small  interval  con¬ 
taining  the  circled  point  on  the  x-axis  into  the  x-axis  is  a  contraction 
mapping  with  respect  to  the  length  of  interval .  We  now  assume  that  the 
actual  system  has  a  transition  such  as  that  from  case  (c)  to  (e)  and 
then  to  case  (d) ,  as  shown  in  Figure  3.  Figure  3(c),  (e)  and  (d)  implies 
that,  as  k  increases,  A2  and  D2  approach  G2  ,  while  B2  and  D2 
approach  each  other  more  slowly.  Then,  the  local  attractor  becomes  narrower 
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and  finally  the  double-loop  local  attractor  case  (e)  appears  and  lasts  until 
case  (d)  occurs.  It  becomes  more  probable  that  the  interval  for  contraction 
mapping  appears  as  dx/dx  decreases  in  the  local  attractor  (indicated  bv 
a  solid  line)  in  Figure  4(a),  which  'anpens  as  the  local  attractor  becomes 
narrower.  Thus,  as  <  increases,  limit  cycles  having  a  small  number  of 
loops  appear.  From  Figure  4(a),  a  single-loop  limit  cycle  must  apnear  after 
a  double-loop  cycle.  It  cannot,  however,  exist  in  the  local  attractor  of 
type  (c) .  This  may  be  the  reason  why  the  interlinked  double  limit -cycles, 
each  of  which  is  a  single- loop,  are  not  observed  in  the  actual  system. 

In  Figure  4(c),  we  also  show  the  graph  of  values  x  of  x  on 
r„  (resp. ,  T_  =  (X  :  z=0  ,  x<  0  ,  y  >  0))  with  respect  to  the  previous  values 
of  x  on  T_  (resp.,  T+)  for  Yj=0.2  and  0.45  .  Figure  4(a)  and  (c) 
show  that  the  return  maps  are  ’folded’  10  for  yt  e  [0.2,0.489].  The  at¬ 
tractor  produced  by  such  a  map  is  called  a  ’chaotic  attractor’.10  In  this 
attractor,  there  may  exist  stable  periodic  solutions.  If  they  exist,  the 
attractor  is  not  a  strange  attractor.  Here,  by  a  strange  attractor,  we 
mean  a  positive  limit  set  of  the  integral  curves  o^  a  differential  equation 
which  is  neither  empty  nor  an  equilibrium  set  nor  a  closed  orbit.  Since 
the  size  of  the  region  of  attraction  for  the  stable  periodic  solution  is 
usually  small,  a  small  perturbation  will  expel  the  trajectories  from  the 
region.  Hence,  even  if  stable  periodic  solutions  exist  in  the  attractor, 
we  actually  cannot  distinguish  a  chaotic  attractor  from  a  strange  attractor. 
We  cannot  prove  here  whether  our  system  has  a  stable  periodic  solution  for 
Y !  «  [0.2,0.489],  Hence,  we  do  not  know  whether  the  attractors  obtained 
in  the  experiment  are  strange  attractors. 

If  the  attractor  A  of  a  system  is  ergodic  and  u(X)  is  a  defined 


measure  on  A  ,  then  for  any  smooth  function  f  , 
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rT 

f  (X)du(X)  =  lin  l  j  dt  £[X(t)]  . 
A  T-«o  h 


(4.4) 


Also,  if  the  time  evolution  of  the  system  is  mixing,  then  for  any  smooth 
functions  f  and  g  , 


«(T)»fg-[  f[X(t)]  g  CX0)du(X.) 

J  A 

=  f  g  -  lim\  fT  dt  f[X(t+r)]  g  [X Ct) ] 
T*»  1  Jo 


(4.5) 


converges  to  zero  as  t  -*■«  .  This  property  means  that  the  time  evolution 
of  the  system  trajectories  is  highly  sensitive  to  the  initial  condition. 
The  solution  with  such  a  property  is  defined  as  a  turbulent  solution.1’5 
We  do  not  know  whether  our  system  has  the  property  of  ergodicity  or  mixing. 
In  order  to  obtain  some  idea  on  the  sensitivity  of  our  system  to  initial 
conditions,  we  have  calculated  the  autocovariance  functions  R(t)  of  x 
for  Yj  =0.2,  0.45  and  0.49  (see  Figure  5).  Although  we  cannot  draw 
any  definite  conclusions  from  the  results,  they  suggest  that  R(t)  conver¬ 
ges  to  zero  as  t-*-°°  and  the  system  has  the  mixing  property. 

We  shall  consider  whether  the  solution  of  system  (3.6),  with  a 
specified  initial  condition  which  is  not  on  the  x,y  or  z-axis,  is  uni¬ 
formly  bounded  on  the  time  interval  [0,°°)  .  We  have  calculated,  for  some 
initial  conditions,  the  ratio  P(x)  of  the  sum  of  the  time  duration  when 
x(t)  exceeds  x  for  te  [0,T)  and  T  ;  namely,  P(x)  =  [length  of 
(t  :  x(t)  *x  5  0  *  t  <T}]/T  .  In  Figure  6,  we  have  plotted  only  those  values 
of  x  for  which  P  is  insensitive  to  the  initial  conditions.  As  T  in¬ 
creases,  the  set  of  such  values  of  x  also  increases  in  size.  Hence,  it  is 
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likely  that,  as  T-*-»  ,  P  approached  ,  which  is  independent  of  the 

initial  conditions.  If  this  is  true,  then  Pm  (x)  t  0  at  finite  values 
of  x  »  for  we  can  choose  an  arbitrarily  large  initial  condition  whose 
corresponding  trajectory  also  has  P^  (y)  •  Thus,  if  P^  exists,  the 
solution  of  (3.6)  with  any  initial  condition  is  not  uniformly  bounded  on 
the  time  interval  [0,°°)  . 

If  P^  exists  and  decays  faster  than  1/y3  ,  then  from  the  physical 
point  of  view,  we  can  state  as  follows:  the  wave  energy  can  be  arbitrarily 
large  over  a  finite  time  interval,  but  its  time  duration  is  so  short  that 

fT  r 

the  time -averaged  energy  lim  (1/T)  x2(t)dt  =  X2pOT(x)dX  is  finite. 

T-*»  o  Jo 

Then,  the  unboundedness  of  system  (3.6)  does  not  contradict  the  physical 
requirement  that  the  energy  of  a  system  must  be  finite. 

We  must  show  that  the  assumptions  made  in  Section  2  in  deriving 
(3.6)  are  consistent  with  the  numerical  results,  or  that  the  assumptions 
are  valid  for  the  observed  values  of  ||X(t)||  .  Although  this  is  true, 

it  is  still  possible  that,  at  some  t  which  is  not  realized  in  the  experi¬ 
ment,  ]|X(t)  ||  is  so  large  that  the  assumptions  are  violated.  If  the 
value  of  ||X(t)  ||  violating  the  assumptions  is  much  larger  than  the  ob¬ 
served  value  of  || X (t)  |{  ,  the  violation  of  the  assumptions  rarely  happens 
and  the  system  (3.6)  is  valid  most  of  the  time. 

Nonlinear  terms  in  (1.1)  produce  only  small  shifts  in  frequency, 
since  |VA2|  or  |VA3|  s  yjxl  <<  ojj  for  the  observed  values  of  |x(t)| 

=  || X (t)  ||  .  The  electrostatic  potential  of  plasma  wave  '  1  *  is  given 


by  <t>j  *  (Y3/V(2k12/3TrUpe)}x  .  Hence,  the  energy  density  per  unit  wave 
number  at  k3  ,  which  is  comparable  to  that  at  k3  ,  is  given  by 
Ej  =  kj 2<>1  2/8tt  =  10’ 18  x2  (erg).  This  is  much  smaller  than  the  kinetic 


energy  of  the  electrons  E^n  =  6.15  *  10' 13 (erg)  for  the  observed  values 
of  x(t).  Hence,  vT  is  constant  since  waves  *1*  and  ’3  are  only  plasma 
waves  excited  in  our  model.  Therefore,  y2  ,  y,  and  V  are  constant  for  the 
observed  values  of  x(t) .  The  bounce  frequency  of  trapped  electrons  is 

i/  i/ 

given  by  Ug  =  (ekj2 |<frx  |/mg)  /2*  3.37xl04|x|  2  (sec-1).  The  beam  particles 
pass  through  the  plasma  in  L/v^  =  1.8  * 10'7  (sec)  ,  which  is  much  smaller 
than  the  trapping  time  l/uig  =  2.97xl0'5|x|  1  (sec)  for  the  observed  values 

of  x(t).  Furthermore,  the  velocity  distribution  function  of  the  beam  is 
time- independent  at  the  surface  of  the  plasma.  Hence,  the  velocity  distri¬ 
bution  function  of  the  beam  can  be  approximated  to  be  independent  of  time 
and  space  in  the  plasma  slab.  Thus,  ys  is  constant  for  the  observed  values 
of  x(t). 


V.  CONCLUDING  REMARKS 

Vie  found  that  the  behavior  of  the  system  of  resonant  three-wave 
interaction  with  linear  growth  and  damping  terms  is  completely  different 
from  that  without  such  terms  or  that  approximated  by  a  two-wave  model. 
Although  it  is  not  proved  mathematically  that  the  system  (3.6)  has  the 
statistical  properties  of  ergodicity  and  mixing,  numerical  results  show 
that  the  solution  (3.6)  is  very  sensitive  to  its  initial  condition.  This 
suggests  that  a  turbulent  state  can  be  produced  directly  by  a  linearly 
unstable  wave  through  the  interaction  with  two  linearly  damped  waves.  We 
note  that,  in  such  a  mechanism,  neither  any  additional  input  of  random¬ 
ness  nor  any  interaction  between  a  large  number  of  waves  is  necessary  to 
produce  a  turbulent  state. 
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System  (3.6)  is  a  highly  simplified  model  with  many  assumptions.  In 
an  actual  plasma,  however,  many  wave-wave  and  wave-particle  interactions 
must  be  considered.  Very  little  is  known  about  the  chaotic  property  of  a 
system  with  dimension  m  >  4  .  It  might  be  completely  different  from  that 
of  a  three-dimensional  system.  Hence,  at  present,  it  is  difficult  to  cor¬ 
relate  the  numerical  results  obtained  for  system  (3.6)  with  the  turbulence 
in  an  actual  plasma. 
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FIGURE  CAPTIONS 


Fig. 1:  Sketch  of  the  surface  S  where  d(y/z)/dt  *  0  ,  the  sets 

Wj  and  W2  separated  by  S  and  the  boundary  3r  of  the 
set  r. 

Fig. 2:  Projections  of  the  phase-space  trajectory  of  equation  (4.1)  onto 

the  (x,z)-  and  (x,y) -planes  for  te  [0,3,276)  and  various 
values  of  Yj  :  (a)  Yi  “ 0*2 »  00  Yi  =  0.45,  (c)  Yi  =  0.489, 

(d)  Yi  =  0.49,  (e)  Yi  =  0.S,  (f)  Yi  =  0.6;  tr  =  transient. 

Fig. 3:  Sketch  of  the  local  attractors  of  the  idealized  system. 

Fig. 4:  Return  mappings  of  the  actual  system. 

Fig. 5:  Autocovariance  functions  R(t)  of  x(t)  obtained  from 

32,768  data  points  for  te  [0,6,554). 

Fig. 6:  The  ratio  P(x)  of  the  sum  of  the  time  durations  when 

x(t)  >  x  for  [0,T)  and  T;  T  *  6,554  . 
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ABSTRACT 


The  stability  of  planar,  mode -converted,  lower-hybrid  solitons 
to  transverse  perturbations  is  investigated.  Two  classes  of  modes  are 
found.  In  the  lowest  order,  the  stability  of  each  class  is  related  to 
a  particular  term  of  the  equation  governing  the  nonlinear  evolution  of 


these  waves. 
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I.  INTRODUCTION 

It  is  well  known  that  one -dimensional  Langmuir  solitons  are 
modulationaily  unstable  to  transverse  perturbations . 1  As  a  result  of 
the  instability,  the  soliton  breaks  up  into  collapsing  bunches  with  longi¬ 
tudinal  and  transverse  dimensions  of  the  same  order.2’8  There  is  much 
interest  in  electrostatic  plasma  waves  with  frequencies  near  the  lower- 
hybrid  resonance  frequency  as  a  possible  method  of  heating  tokamak  plasmas 
According  to  linear  theory,  cold-plasma  waves,  excited  by  an  rf  source 
at  the  surface  of  a  magnetically  confined  plasma,  propagate  in  resonance 
cones3  into  the  interior  of  the  plasma  until  they  reach  the  lower-hybrid 
resonance  layer,  where  they  may  mode-convert  to  a  slow  hot-plasma  wave.4,s 
Nonlinear  ponderomotive  effects  on  the  cold-plasma  wave6 ,  as  well  as  the 
mode -converted  lower-hybrid  wave7,  have  been  examined  in  a  two-dimensional 
geometry.  Recently,  the  effect  of  a  third  dimension  on  the  cold -plasma 
wave  was  studied  and  it  was  found  that  perturbations  transverse  to  the 
soliton  structure  were  unstable . 8  An  equation  governing  the  nonlinear 
dynamics  of  mode -converted  lower-hybrid  waves  in  three  dimensions  has  also 
been  presented. 8 

In  this  paper,  we  address  ourselves  to  the  stability  of  a  mode- 
converted,  lower-hybrid,  planar  soliton  to  perturbations  in  the  transverse 
direction.  We  find  that  these • solitons  are  also  unstable. 


2 


II.  EXISTENCE  OF  UNSTABLE  MODES 

The  equation  governing  the  nonlinear  evolution  of  mode-converted, 
lower-hybrid  waves  is,  according  to  Ref.  [9], 


1  Jr,j;  +  ao'^2  *  +  6o  lW4'  +  Yo'^T  ^  +  5o 


Jy  14*  |  2|^  *  0 


(1) 


where 


a  -  k.) 
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K  5  1  +  K)  ^  "  K>  *  K||  *  1  -  Ca>  AO2 

J‘  *e  ®  Fi  " 

A0  -  f  CV/*)2Cv./U)2  +  j  Cw  /«c  )2(v>c  )2 

Fi  pe  e  e 

L2  u) 

Pi 


CO  <0 

P  C 
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(U2-U  2)  <o(u>2-!U  2) 

c.  ci 


in  which  to  and  to  denote  the  ion  and  electron  plasma  frequencies, 
pi  pe 

~c  and  toc  are  the  electron  and  ion  cyclotron  frequencies,  and 
e  i 

v^2  *  2T./m^  ,  vg2  *  2Tg/me  .  In  Eq.(l),  ip  is  the  normalized  field 

i/ 

amplitude  [XJ/4no(T..  +  T  )Aq]  ri  $  associated  with  the  electrostatic  poten¬ 
tial  i  of  'the  wave  .which  has  the  form 


$  *  4>  expfiCKj/A^  ^  X  -  itot]  . 


(2) 


/ 
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(3) 


W 


i  jf'b  *  +  <|^|21J/  *  ct(~  |^|2)  ^  +  6^5-i(>  *0  ,  (5) 

which  is  the  usual  nonlinear  Schrodinger  equation  with  two  additional  terms. 
For  a3  8*0  ,  Eq.  (5)  is  known  to  have  the  solution 

^o  3,n/k'A  secil CAC) exp [ - x<aJqt  +  ikQz] 

=  <f>(C)exp[-ia>0T  +  ikQz]  (6) 


where 


5  -  Z  -  Vt  ,  V  *  2k  ,  A2  -  k  2  -  u. 

*  O  0  0 


which  is  stable  to  perturbations  in  the  z -direction.  . 

Next,  we  linearize  Eq.(5),  setting  ij>  *  i|>Q  +  6ip  and  assume  that 


-iw  r + ik z 

<5^  *  [F(C,T)cos(kyy)  +  GCC,x)sinCkyy)]e  .  (8) 

Eq.  (5)  then  reduces  to  a  linear  combination  of  terms  proportional  to 
sin(kyV)  and  cos (k^y)  .  Setting  the  coefficients  of  each  of  these 
equal  to  zero  yields 
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The  only  eigenfuction  of  Lq  that  vanishes  at  infinity  has  an  eigenvalue 
of  zero.  Therefore,  according  to  Eq.(lS),  ky =* 0  is  the  only  marginally 
stable  state  corresponding  to  n =  0  . 

Eqs. (15-16)  can  be  reduced  to  the  form 


(LQ  +  Sky2)  [(1^  +  Sky.2 ) E  -  2akyJ2o]n  =  fl2n 
[  (Lx  ♦  Sk  2)E  -  2aq2a](Lo  +  8ky2)?  -  a2c 

For  mode-converted  lower-hybrid  waves,  K  -  1  ,  K^~  e 

where  e2=m/m.  .  This  implies  that 

Q  X 


akyfrVA2  «  1 


(17a) 

(17b) 
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(18a) 

(18b) 


aky$2/ic!!)2  «  1  ,  (18c) 

where  k_  is  the  component  of  the  propagation  vector  in  the  original  coor¬ 
dinate  system.  Also,  k,2/k  2  «  e2  ,  and  ky/k2  £  1  . 7,9  For  values 

of  k  in  this  range,  Eq.  (17a)  can  be  solved  by  expanding  n  and  ft2  in 
the  form  of  a  perturbation  series. 


q  *  riQ  +  +•  *  ■  ,  ft2  a  +  +'  ’ '  »  ©be. 


In  lowest  order,  Eqs. (17a, 17b)  reduce  to 

L0Lln0  ’  0  •  Wo  '  0  •  tl9) 
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In  order  to  determine  n  and  t  ,  consider  the  functions,1 

oo  ’ 


+  9.  ,  3  ,  + 

uo  "  '  5a5"*  ’  u0  ■  2  at *  ’  vo  3  -  v0 


-£4> 


(20) 


According  to  Eq.(14), 


L  v  =  0  • 

o  o 


(21) 


3y  differentiating  Eq.  (14)  with  respect  to  £  and  A2  ,  we  find  that 


L  v  =  u 
oo  o 


.  +  + 
L,u  »  v 
I  0  0 


Vo  ■  0 


(22) 


Eqs.  (21-22)  imply  that 


LL.U" 
oi  o 


-  hVo  ■" 


The  only  solutions  of  Eqs.  (19)  that  vanish  at  infinity  axe  u*,v*  . 


(23) 

Therefore , 


c,u  +  c,u 
1  o  2o 


?o  *  Cl'v0  +  c2,vo 


are  the  proper  eigenfunctions  of  Eqs. (19),  where  c^,  c2>  c^' ,  c^'  are 

constant  vectors.  The  parity  of  u+,v*  is  even,  while  that  of  u’,v  is 

0  0  0  0 

odd  with  respect  to  £  . 

In  first  order,  Eq. (19a)  takes  the  form 


L0Llril  +  -  L02^2?In0  3  Exn0  . 

.'e  define  the  scalar  product  of  two  functions  as  follows: 


(24) 


<UIV>  2 


d£uv 


J  -co 


Taking  scalar  products  of  Eq. (24)  with  v~  yields 
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2  +,  + 


lVTolVcl  ■  «I<voluo>cl 


(25) 


and 


[gkJ<uo|uo>E  -  2aky<u^62!u‘>a]c2  =  Ci<v^|u^>c2  ,  (26) 

where  we  have  used  Eqs.  (21-23)  and  the  self-adjoint  property  of  L  and  L, . 


We  note  that 
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Therefore, 


3+2  S  -  -  48k  2A2  . 


(28) 


The  eigenvalue  ^  is  obtained  from  Eq.  (26)  by  solving  the  detemin- 
antal  equation 


sy<“;iv-.q<v;iv 

-  2aky<u^|<J)2|u‘  > 

2aI^<uol^2l%> 

®V2<uoK>-  5l<voluo> 

(29) 


rcm  which  we  find 
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hmT*  VA2^i§  < 


(30a) 


-  ±(a2  +  b2)  e 


%  .±ia/2 


(30b) 
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where 

a  =  i  3ky2A2  ,  b  =  Ctk^  ^  ,  3  *  arctan(b/a) 

with  growth  rates 

Y,  3  2  k  A*/6 

'  *+•  y 

16  "V4 

Y-  =  IT  — 

2  2 

Since  the  mode-conversion  region  is  excluded  (k„/kx)|Kj|  |  «  1  ,  from  which 
we  find  that  3kv</ctA2  f  e2  ,  so  that 

Y_  =  VlT  • 


III.  CONCLUDING  REMARKS 

According  to  Eqs. (28,30),  unstable  modes  of  either  parity  exist.  In 
first  order,  positive  parity  states  are  driven  unstable  by  the  term 
Y0(32/3Y2) .  If  6-  were  absent,  Eq.(l)  would  be  of  the  same  form  as 
the  equation  governing  the  nonlinear  evolution  of  Langmuir  waves  with  ion 
inertia  neglected.  It  is  this  term  which  is  responsible  for  the  instability 
of  Langmuir  solitons.  Instabilities  associated  with  negative  parity  modes 
are  due  in  lowest  order  solely  to  the  additional  nonlinear  term 
^[(s/aY)  '-yj2]  y  ,  which  affects  the  stability  of  positive  parity  states  in 
second  order.  The  relative  effects  of  these  two  terms  can  be  judged  by 
comparing  growth  rates  of  positive  and  negative  parity  modes.  We  find  that 
yj'(_  $  10  e  < .  1  in  regions  which  are  not  mode-converting,  which  indicates 
that  the  term  5  [(3/3Y)  |u(2]  y  is  strongly  destabilizing.  We  therefore 
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expect  large  local  values  of  the  electric  field  to  exist  as  the  wave  prop¬ 
agates  away  from  the  lower-hybrid  layer,  thus  causing  the  soliton  to  break 
up  into  many  bunches,  which  move  apart,  spreading  energy  throughout  the 
plasma. 
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ABSTRACT 


Chaotic  oscillations  or  strange  attractors  in  bilinear  systems 
with  linear  feedback  controls  are  explored,  A  simple  sufficient  condi¬ 
tion  for  nonexistence  of  such  oscillations  is  given.  Also,  the  con¬ 
struction  of  ellipsoidal  bounds  for  the  strange  attractor  is  discussed. 
The  results  are  applicable  to  the  well-known  Lorenz  equation  and  other 
eq'iations  which  are  known  to  have  chaotic  oscillations. 
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I.  INTRODUCTION 

Recently,  it  was  discovered  that  the  nonlinear  mathematical  models 
of  a  number  of  physical  processes  exhibit  chaotic  nonperiodic  oscillations. 
A  simple  example  is  the  Lorenz  model  [1]  for  thermal  convection  in  a  fluid 
layer  given  by 

xil  °(x2  “  xi^ 

x2  <*  rXj  -  x2  -  xxx3  ,  (1) 

x3  -bx3  +  XjX2 

where  o,r  and  b  are  constant  parameters.  It  was  found  that  for  certain 
ranges  of  these  parameters,  the  trajectories  of  (1)  in  the  state  space  are 
attracted  to  a  nonempty  set  which  is  neither  an  equilibrium  set  nor  a 
periodic  orbit.  On  this  set,  the  trajectories  exhibit  chaotic  nonperiodic 
oscillations.  Ruelle  and  Takens  [2]  called  such  a  set  a  "strange  attrac- 
to-".  Oscillations  having  similar  chaotic  behavior  were  discovered  in 
other  models  such  as  those  for  chemical  turbulence  [3],  population  dynamics 
[4],  disk  dynamos  [5]  and  plasma  turbulence  [6],  A  remarkable  feature  of 
such  oscillations  is  that  their  apparent  chaotic  or  random  behavior  is 
intrinsic  to  the  model  which  is  completely  deterministic  in  nature. 

In  this  paper,  we  consider  a  bilinear  system  of  the  form: 

dx/dt  ■  Ax  +  uBx  (2a) 

defined  on  the  n-dimensional  Euclidean  space  En  whose  inner  product  is 
denoted  by  (*,•).  A  and  B  are  real  constant  n*n  matrices,  and  u  is  a 
linear  feedback  control  given  by 

u  -  -K^c,x),  (2b) 

where  c  is  a  nonzero  vector  in  En,  and  K  is  a  real  number  corresponding 
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to  the  feedback  gain.  We  note  that  the  Lorenz  equation  (1)  can  be  re¬ 
written  in  the  form  of  (2)  with 

■  -a  a  01  TOO  O' 

A-  r  -1  0  ;  B  -  0  0-1  ;  c  -  (1,0, 0)T  (3) 

.0  0-bJ  L010. 

and  K  *  -1. 

Although  extensive  studies  have  been  made  on  bilinear  control  systems 
[7]-[10],  most  of  the  attention  has  been  focused  on  optimal  control,  vari¬ 
ous  system  properties  such  as  controllability  and  observability,  and  reali¬ 
zation  theory.  Very  little  work  has  been  devoted  to  the  self-oscillation 
of  bilinear  systems  with  linear  or  nonlinear  feedback  controls.  Here,  we 
explore  the  possibility  of  existence  of  chaotic  nonperiodic  oscillations 
in  system  (2) .  We  begin  by  studying  the  stability  of  the  equilibrium 
states  of  (2),  and  obtaining  some  estimates  for  the  solutions.  Then,  the 
existence  or  nonexistence  of  chaotic  oscillations  is  explored. 

II.  STABILITY  OF  EQUILIBRIUM 

Assuming  that  A  is  nonsingular,  the  equilibrium  set  E  of  system  (2) 
consists  of  all  fixed  points  of  the  mapping  H  defined  by 

H(x)  -  K  ^c,x)  A  1 Bx.  (4) 

Obviously,  the  zero  vector  belongs  to  E.  Let  xg  be  a  nonzero  vector  in 
E.  Then, 

|xj«  |K||<c.*>||A-1B«J  *  |K|||cj|||A-lB||||xeP.  (5) 

Hence,  any  nonzero  equilibrium  state  must  lie  outside  or  on  the  sphere 
with  radius  p  -  (|k|Jc|||a  1b||)  l.  Also,  since  H  is  an  even  operator 
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(i.e.  H(-x)  *  H(x)  for  all  x£En),  hence  for  any  nonzero  xg  G  E,  ~xe$  E, 
Now,  consider  Che  equation 

x  -  H(x)  (6) 

for  determining  the  equilibrium  states.  Since  K(c,x)  is  a  real  number, 
any  nonzero  solution  x&  of  (6)  must  be  an  eigenvector  of  A  *B  correspond¬ 
ing  to  some  nonzero  real  eigenvalue.  Moreover,  ^c,x^  ^0.  Hence,  any 
nonzero  equilibrium  state  xg  must  belong  to  n(A  *B)  f|{x6En:  (c,x^  /  0}, 
where  n(A  lB)  denotes  the  null  space  of  A  1B.  Let  v^  be  an  eigenvector 
of  A  iB  associated  with  the  nonzero  real  eigenvalue  A^.  Then,  its  cor¬ 
responding  equilibrium  state  xg  has  the  form  ctv^,  where  the  scalar  a  can 
be  determined  by  substituting  av^  into  (6),  i.e. 

avi  “  K(c»avi)  A_IB(avt)  »  AiKa2(c,vi>vi 
or 

a  -  l/(XiK  (c.v^),  (7) 

provided  that  K  (c.v^)}*  0. 

Evidently,  if  A  lB  has  n  distinct  nonzero  real  eigenvalues  A^  with 
corresponding  eigenvectors  v^  then  the  equilibrium  set  E  is  given  by 

E  *  {0}  U  (vi/(AiK  (c.v^)):  K(c,vi)  i  0,  i-l,...,n}  .  (8) 

Since  {v^,i«l, . . .  ,n}  is  linearly  independent  and  c  +  0,  there  exists  at 
least  one  i,  say  i-k,  such  that  ^c,v^  +  0.  Thus,  for  any  nonzero  K,  E 
has  at  least  two  and  at  most  (n+1)  distinct  points.  In  the  case  where 
A  lB  is  a  simple  linear  transformation  with  repeated  real  eigenvalues, 
say  A^^  ■  A  with  eigenvectors  v^,  i€  I  €  {1, . . .  ,n} .  Then,  for  K  /  0  and 
c^  [span{v^,i6 I}]-*-,  E  has  an  Infinite  number  of  equilibrium  points, 
since  any  nontrivial  linear  combination  v  of  such  v^'s  is  again  an  eigen- 
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vector  of  A  lB  with  (c,v)  i  0. 

To  study  the  stability  of  an  isolated  equilibrium  state  xg,  we  con¬ 
sider  the  following  linearized  system  about  xg: 

d6x/dt  “  J^(xe>6x,  (9) 

where  Jf(xg)  is  the  Jacobian  matrix  of  f(x)  -  Ax  -  K(c,x)Bx  at  xg  given 
by 

Jc(x  )  -  A  -  KB{(c,x  >  I  +  x  cT>  .  (10) 

f  e  e'  e 

Obviously,  the  stability  of  the  origin  is  determined  only  by  the  eigen¬ 
values  of  A.  At  a  nonzero  equilibrium  state  xg  ■  v^/CA^K  ^c,Vj>),  we 
have  A  *Bx  ■  Ax  or  Bx  ■  Av  /  (K  ^c,v^).  Thus,  J,(x  )  can  be  re- 
written  as 

Jf(xe)  -  A  {I  -  <c,vi>-1vicT}  -  A“*B  (11) 

which  is  independent  of  the  feedback  gain  K.  This  is  evident  from  the 
fact  that  for  Kf*0,  we  may  introduce  the  scaling  x=x/K.  Thus, (2)  is  equi¬ 
valent  to  the  system 

dx/dt  -  Ax  -  ^c,x^Bx.  (2') 

III.  ESTIMATES  FOR  SOLUTIONS; INVARIANT  SETS 

First,  we  derive  a  few  elementary  estimates  for  the  solutions  of  (2) 
under  various  assumptions  on  A  and  B. 


(B  +  B^)  respectively,  and  h  •  max{M.},  | A |  -  max{|A  |).  Assume  that 

A  i  j  J 

v  -  | K 1 1 c  || | X |  >  0.  Then,  a  solution  of  (2)  starting  with  x  at  t  ■  0 

satisfies  one  of  the  following  estimates: 


(i)  for  P“0: 

1  x(t)  ||  <  ||  xq  1|(1  -  v  ||xqJ  t/2)  1  for  all  t  6[0,2(|xqJv)  j[; 

(ii)  for  y>0: 

|x(t)||  <  II  xq  ||  [  <u  +  v|xj  )exp(-yt/2)  -vjxjr1 
for  all  where 

■  2y-1£n[l  +  u(v|jxQ  jj)  1  ] ; 

(ill)  for  U<0  and  f |]  £  jyj/v,  |x(t)l|  ^  1xq||  for  all  t^O.  Moreover, 
||x(t)  ||  -*■  0  as  t-*-°°. 

Proof:  Let  V(t)  -  ||x(c)||2.  By  direct  computation 

dV/dt  -  <x,(A  +  AT)x>  -  K<c,X><x,(B  +  BT)x> 

<  yV  +  |K||c|||x|||<x,(B  +  BT)x>|  <  (y  +  VV2)V  -  h(V)  (12) 

with  V(0)  -  SXQI|2* 

For  y>0,  h  is  a  strictly  monotone  increasing  function  of  V. 

Consequently,  V(t)  -  ||x(t)||2  4  w(t),  where  w  is  the  solution  of 

1 

dw/dt  ■  (y  +vw2)w,  w(0)  ■  llx0ll  •  (13) 

A  A 

Solving  the  above  equation  for  each  case  with  y»0  or  p>0  leads  directly  to 
the  estimates  (i)  and  (ii)  .  When  p<0,  dV/dt<0  for  0<V  <  ( |  y  |  /  V) 2 .  Re¬ 
sult  (ii)  follows.|| 

T 

In  the  special  case  where  B  is  skew-symmetric  (i.e.  B--B  ),  we  have 
the  exponential  estimate: 

||x(t)||  <  | xq|| exp (y t/2)  for  all  t*0,  (14) 

which  implies  that  no  solution  has  finite  escape  time.  Finally,  when  both 
A  and  B  are  skew-symmetric,  |x(t)|II|x0l  for  all  t  or  the  solution  remains 
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on  the  sphere  with  radius  ||  xq ||  at  all  time. 

Now,  we  give  lower  bounds  for  the  solutions. 

y  T 

Proposition  2:  Let  y  denote  the  minimum  eigenvalue  of  (A  +  A  )  and 
|X|,v  be  as  defined  in  Proposition  1.  Then  the  solution  of  (2)  starting 
from  xq  at  t»0  satisfies  one  of  the  following  lower  bounds: 

(i)  for  y-0, 

||x(t)||  >  |xo||/(l  +  v|xJt/2)  for  all  t;>0;  (15) 

(ii)  for  yj*0, 

lx(t)||  >  m||xo||/[(u  -  V  ||xo|)exp(-yt/2)  +  v|xj|]  for  all  t>0.  (16) 

Proof:  Let  V(t)  -  ||x(t)|2  Then, 

dV/dt  >  yV  -  K<c,x)<x,(B  +  BT)x)  >  (y  -W2)V  (17) 

with  V(0)  -  || ||2.  Consequently,  V(t)  »  ||x(t)||2  »  w(t),  where  w  is  the 
solution  of 

1 

dw/dt  ■  (y  —  V w2 ) w,  w(0)  *  IX0I|Z*  (18) 

Solving  (16)  leads  to  estimates  (15)  and  (16). || 

We  observe  that  for  y>0  and  y<0,  the  lower  bound  (16)  tends  to  y/v 
and  zero  respectively  as  t  ■+■».  Also,  the  lower  bound  (15)  tends  to  zero 
as  t  -*•  <“. 

Now,  we  consider  the  existence  of  linear  subspaces  in  En  which  are 
invariant  sets  of  the  system. 

Proposition  3:  Let  A  represent  a  simple  linear  transformation  op  En 
with  distinct  real  eigenvalues  X^  and  their  corresponding  eigenvectors  v^, 
i“l . .  Let  JC{l,...,n}  denote  the  index  set  such  that  (X^iGJ)  re¬ 

presents  all  the  uncontrollable  modes  of  (A,B)  and  all  the  uncbservable 
modes  of  (A,c  ).  Then,  span{v^;i€J}  is  an  invariant  set  of  (2). 
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-1 


Proof:  Let  T  =  [v, ].  Since  A  is  simple,  T  ‘exists  and  T  *AT 
-  1 '  '  n 

A  =  diag[A^,...,A  ].  Let  x  =  Tz.  Then,  (2)  is  transformed  into  the 
following  form: 

dz/dt  -  Az  -  K<c,Tz>  T_1BTz. 


(19) 


If  A^  is  an  uncontrollable  mode  of  (A,B),  then  the  j-th  row  of  T  B  must 

/V  p* 

be  zero.  Also,  if  A^  is  an  unobservable  mode  of  (A,c  ),  then  the  j-th 
T 

column  of  c  T  must  be  zero.  In  either  case,  the  equation  corresponding  to 
the  j-th  mode  reduces  to  dz^/dt  «  ^jzj*  TLUS»  the  subspace  spanned  by 
~  is  an  invariant  set  of  (2),  and  the  desired  result  follows. || 

For  a  general  linear  transformation  A,  we  can  introduce  the  usual 
canonical  decomposition  of  the  state  ipace  of  the  linear  system:  dx/dt  = 

Ax  +  Bu,  y  =^c,x^  according  to  its  controllable,  uncontrollable,  observa¬ 
ble  and  unobservable  modes.  It  is  easy  to  see  that  the  subspace  cor¬ 
responding  to  the  uncontrollable  and/or  unobservable  modes  of  the  fore¬ 
going  linear  system  Is  an  invariant  set  of  (2). 

IV.  CHAOTIC  OSCILLATIONS 

The  existence  of  chaotic  oscillations  or  strange  attractors  depends 
on  the  manner  in  which  the  stable  and  unstable  manifolds  intersect  with 
each  other.  Given  a  system  such  as  (2),  suppose  we  could  find  all  the 


bounded  invariant  manifolds  of  the  system,  then  we  could  seek  chaotic  oscil¬ 
lations  by  deleting  those  invariant  manifolds  which  correspond  to  the 
equilibrium  set  and  periodic  orbits.  The 

remaining  ones,  if  they  exist,  may  consist  of  solutions  which  are  almost 
per iodic, peeudo-random  or  chaotic  functions  of  time.  The  distinction 
between  the  almost  periodic  solutions  from  the  pseudo-random  solutions 
can  be  accomplished  by  investigating  the  asymptotic  properties  of  their 
correlation  functions,  provided  that  these  solutions  are  known.  The 


foregoing  approach,  although  conceptually  simple,  represents  a  formidable 
task  from  the  analytical  and/or  computational  standpoints.  So  far, cha¬ 
otic  oscillations  in  certain  nonlinear  systems  were  discovered  through 
numerical  computation  and  bifurcation  analysis.  It  is  desirable  to  have 
readily  verifiable  analytical  conditions  for  determining  the  existence  of 
chaotic  oscillations  directly  in  terms  of  the  right-hand- sides  of  the  sys¬ 
tem  equations.  Lacking  such  results  at  the  present  time,  we  shall  re¬ 
strict  ourselves  to  the  less  ambitious  tasks  of  establishing  simple  suffi¬ 
cient  conditions  for  the  nonexistence  of  chaotic  oscillations,  and  obtain¬ 
ing  bounds  for  the  amplitudes  of  the  chaotic  oscillations  when  they  exist. 

T 

Theorem  1:  If  A+A  is  negative  definite}  B  is  skew-symmetric,  and 
c  6  n(A  1B),  then  system  (2)  has  no  periodic,  almost  periodic  or  chaotic 
oscillations. 

Proof:  Since  any  nonzero  equilibrium  state  of  (2)  must  belong  to  r)(A-1B)"L 
n{x6En:(c,x)  i*  0},  hence  the  zero  state  is  the  only  equilibrium  state  if  c 
1(A~‘  B) .  Since  A+A  is  negative  definite,  <x,  (A+AT)x>$)J ||x||2for  all  x£E  .where 
is  the  minimum  eigenvalue  of  A+A1" .  Since  B  is  skew-symmetric,  x(t)  has  a  decaying 
exponential  upper  bound  giver,  by  (14) ,  implying  that  the  origin  is  asymp¬ 
totically  stable  in  the  large.  Hence,  no  periodic, almost  periodic  or 
chaotic  oscillations  can  exist. | | 

Assuming  the  existence  of  chaotic  oscillations,  we  proceed  to  con¬ 
struct  an  ellipsiodal  domain  fl  in  En  which  contains  the  invariant  mani¬ 
fold  generated  by  the  chaotic  oscillations.  This  invariant  manifold 
does  not  contain  any  equilibrium  states. 

Consider  a  quadratic  form  in  x  given  by 

V(x)  -  <x-x,Q(x-x)>,  (20) 

^  II 

where  x  is  a  constant  vector  in  E  and  Q  is  a  positive  definite  symmetric 


matrix.  We  shall  make  use  of  V  to  establish  conditions  under  which  ft  is 
an  attractor  of  (2)  in  the  sense  chat  every  trajectory  initiated  from 
the  exterior  of  ft  eventually  enters  ft  at  some  finite  time  t^>0  and  remains 
in  ft  for  all  t>t^,or  tends  to  £1  as  c**.  Clearly, such  a  ft  contains  all  the 
trajectory  points  of  the  chaotic  oscillations  when  they  exist. 

Theorem  2:  Suppose  that  A  is  nonsingular  and  there  exist  a  nonzero 
vector  x£En  and  a  positive  definite  symmetric  matrix  Q  such  that 

(i)  QB  is  skew-symmetric; 

(ii)  G-A^Q-KJA+KP  is  negative  definite, where  P  •  cxTQR+BTQxcT, 
then  there  exists  an  ellipsoidal  set  Q**{x6Etl:l/(x)  ^  a}  which  is  an  at¬ 
tractor  of  system  (2)  containing  the  manifold  generated  by  the  chaotic 
oscillations  when  they  exist. 

Proof :  Consider  dl//dt  given  by 

dV/dt  ■  <x»  (ATQ+<JA)x>  -  K^c,x)(x,  (BT$+<}B)x^  -  2^x,ATQx^ 

+  k(c,x)((x,BTQx>  +  (x.QBx)).  (21) 

In  view  of  condition  (i),  the  second  tern  in  the  right-hand-side  of 
(21)  vanishes.  Rewriting  the  last  term  in  (21)  as  K^x,PxX  (21)  reduces 
to 

dV/dt  »  {x,Gx^  -  2(x,A^QxX  (22) 

From  (ii),  G  1  exists.  We  can  rewrite  (22)  as 

dV/dt  *  (x-xg,G(x-xg)}  -  C,  (23) 

where  xg  -  G  A  Qx  and  C  ■  ^xg,Gxg).  Since  A  is  nonsingular,  Q  is  posi¬ 
tive  definite  and  X  +  0,  hence  xg  0.  Evidently,  at  any  point  exterior 

W  A 

to  the  ellipsoidal  set  ft  ■■  {x&En:  -  (x-x  ,G(x-x  )}<|c|)  ,  dV/dt  <  0. 

s  s 

Also,  since  ft  is  bounded,  there  exists  a  real  number  a  such  that  the  ellip¬ 
soidal  set  ft  *{x£En:  l/(x)  ^  a}  contains  ft.  Consequently,  we  have  dl//dt 
<0  along  any  trajectory  exterior  to  ft,  which  implies  that  ft  is  an  attractor 
of  (2). | | 
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Remarks : 


(R-l)  Condition  (i)  requires  the  existence  ef  a  positive  definite 
symmetric  matrix  Q  such  that 

BTQ  +  QB  -  0.  (24) 


Rewriting  (24)  in  the  usual  form  of  a  linear  equation  Sq  *  0  with 
S  *  BT®I  +  I  ®BT  and  q  «  (q^, . . .  ,q^)^,  where  is  the  i-th  column  of 
Q  and  ®  denotes  the  Kronecker  product,  we  see  that  for  the  existence  of 
a  nonzero  q  or  Q  satisfying  (24),  S  must  be  singular.  Since  the  eigen¬ 
values  of  S  have  the  form  +  Aj  ,  where  A^  and  A^.  are  eigenvalues  of  B, 
there  must  exist  X^  and  X^  such  that  X^  +  X^  =*  0.  This  is  possible  if 
and  only  if  at  least  one  of  the  following  conditions  holds:  (i)  B  is 
singular;  (ii)  B  has  real  eigenvalues  symmetric  about  the  origin,  and 
(iii)  B  has  one  or  more  complex  conjugate  pairs  of  pure  imaginary  eigen¬ 
values.  For  the  special  case  where  B  is  skew-symmetric,  all  its  eigen¬ 
values  lie  on  the  imaginary  axis.  Hence  (24)  has  nontrivial  solutions 
Q.  The  requirement  that  Q  be  positive  definite  imposes  further  restric¬ 
tions  on  B.  Results  for  a  special  form  of  B  were  given  in  [11]  and  [12], 
(R-2)  Condition  (ii)  requires  the  negative  definiteness  of  G.  We 
note  that  if  B  is  singular,  then  for  all  x€n(QB)  (the  null  space  of  QB) , 

<x,Px^  *  0  regardless  of  the  choice  of  x.  Therefore,  in  this  case,  a 

T  T 

necessary  condition  for  G  to  be  negative  definite  is  that  II  (A  Q  +  QA)II 
be  negative  definite,  where  the  matrix  II  represents  the  projection  onto 


n (QB) . 

(R-3)  Given  a  pair  (x,Q)  satisfying  conditions  (i)  and  (ii)  of 

A# 

Theorem  1,  the  smallest  ellipsoidal  set  0  containing  can  be  found  by 
solving  the  following  standard  constrained  optimization  problem:  Maximize 


V  (x)  =  ^x-x,Q(x-x)^  over  the  ellipsoid  3(2*  {x£Ee:<x-x  ,G(x-x  -  C} . 

Although  this  problem  can  be  solved  exactly,  it  is  useful  to  construct  an 

ellipsoidal  set  fl*  (not  necessarily  the  smallest  one)  which  contains  (2. 

L  „  1  i  i 

Let  w  -  Q2 (x-x)  so  that  Vm  |wjj2  and  (2  «{w6En:-(w-Q2  (x  -x)  ,  Q~?GQ~2  (w- 
1  S 
Q2(xg-x))<  }  C |  }  .  Since 

-  1  _i  i  ± 

y||w-Q2(x  -x)||2  $  -<w-QT(x  -x),Q'2GQ~r(w-Q2(x  -x))>  <  |c|,  (25) 

3  S  S 

where  y  is  the  minimum  eigenvalue  of  the  positive  definite  symmetric  matrix 
l  i 

-Q  2GQ  2 ,  we  have 

HI  “  |Q2(x3-x)!i  4  |  w-Q2  (xg-x)  [j  4  (|c|/y)2  (26) 

or 

i  i  A 

M|  4  (|C|/Y)2  +  ||Q2(xs-x)l  £  /a’.  (27) 

Thus,  ihe  ellipsoidal  set 

fl'  -  {x€  !n:  <x-x,Q(x-x)>  ^  a'}  (28) 

contains  (2. 

An  alternate  approach  Is  to  find  a  pair  (x,Q)  satisfying  conditions 
(i)  and  (ii)  such  that  the  ellipsoidal  set  containing  (2  has  the  smallest 
volume.  This  problem  is  not  so  straightforward. 

(R-4)  If  system  (2)  has  invariant  sets  of  the  form  spanCv^  ;i€:J} 
as  described  in  Proposition  3,  then  their  intersection  with  the  invariant 
manifold  generated  by  the  chaotic  oscillations  must  be  empty.  Moreover, 
the  invariant  sets  span{v^}  associated  with  the  positive  eigenvalues 
iGJ  can  have  at  most  the  zero  vector  in  common  with  the  ellipsoidal  set 
(2  in  Theorem  2. 


j 
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V.  EXAMPLES 

First, we  consider  an  equivalent  version  of  the  Lorenz  equation  in  the 
form  of  (2)  with  variable  feedback  gain  K,  and  A,B,c  given  by  (3). 

Assuming  that  b^O  and  a(l-r)f*0,  A  is  invertible.  The  eigenvalues  of  A  lB 
given  by 


A-1B 


(1-r) 

(1-r) 


-n 


-i 


0  -b-1  0 


(29) 


are 


-1 

0,  >-2  3’  2  • 


Thus,  the  equilibrium  set  E  consists  of  the  zero  state  and 

v2/(A2K{c,v2))  *  K  1 (/b(r-l) ,/b(r-l) ,l-r)T, 
v3/(A3K(c,v3>)  -  K  1  (-/b(r-lj  ,->^b(r-l)  ,l-r)T, 

-l 


(30) 


(31) 


where  v2  and  v3  are  eigenvectors  of  A  B  corresponding  to  Az  and  A3  re¬ 
spectively.  Also,  the  eigenvalues  of  A  are  given  by 


\  -[-(l-KJ)+A2]/2,  X2=[-(l+a)-A^]/2,  A3— b 

and  their  corresponding  eigenvectors  are 

ir  f 

2a/ (a-l-A  ) 


(32) 


T  -  [vj.Vj.VjJ 


2a/(a-l+A  ) 
1 


1 

0 


0 

0 

1 


(33) 


,-i 


where  A«(l+a) -4a(l-r) .  If  a  and  A  are  nonzero,  then  T  exists  and 


T-IB  - 


0 

0 

1 


-4arA*“  [a-l-A*]-1 

.L  1  . 

4arA  2  [o-l+A2] 


(34) 


A* 


C  T  -  [2a/(a-l+A2) ,2a/(a-i-A2) ,0] . 
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Evidently,  (A,B)  is  completely  controllable  and  X3  is  the  only  unobserva¬ 
ble  mode  of  (A,cT).  From  Proposition  3,  span{v3}  is  an  invariant  set  of 

the  system.  Since  B  is  skew-symmetric,  the  solutions  have  an  exponential 

a  ,  ,  1 

upper  bound  (14)  with  max{-2b,-(l+o)±[  (1-c)  +(c+r)  ]2}. 

Now,  we  apply  Theorem  2  to  construct  an  ellipsoidal  set  ft  which  con¬ 
tains  the  orbits  of  the  chaotic  oscillations  when  they  exist.  To  simplify 


the  computations,  we  make  use  of  a  function  V  of  the  form  (20)  with  a  di¬ 
agonal  matrix  Q  -  diagUj  t  ,q2  2  ,q , ,  ] .  To  satisfy  condition  (i)  of  Theorem 

**  ~  ~  ~  T 

2,  we  must  have  q22«qj,.  Let  x  -  (Xj.Xj.Xj)  .  Then,  the  matrix  G  in 
condition  (ii)  of  Theorem  2  is  given  explicitly  by 


-2oq,,  oq,  I  +  (r-HCxJ)qJJ 

oq11+(t+^3)q3 ,  -2q}J 


-Kq3Jx2 

0 


L-KqS3*2  0  -Zbc«3  3  J 

We  must  choose  qu,q33>0  and  x  such  that  G  is  negative  definite, 
pedience,  we  set 


(35) 


For  ex- 


*2  “  °.»  *3  "  rq33>/^Kq33)  (36) 

so  that  G  reduces  to  a  diagonal  matrix.  For  O,b>0,  G  is  negative  definite 

for  any  positive  qu  and  q33.  Now,  we  have  dV/dt <  0  at  any  point  exterior 

the  ellipsoidrlSCtft  »  {x€rE*:  -^x~x  ,G(x-x  )>  4  |  C  j  } ,  where 

s  s 

xs  -  G_1ATQx  -  (xl/2,-Oq11x1/(2q33),-(oq11+rq33)/(2Kq33))T,  ^ 

C  -<xs,Gxs>=-  -(l/2)q73{aqllxl(qJ1+Oq11)+b(aq1J+rq33)  /K  }. 

The  remaining  task,  is  to  construct  an  ellipsoidal  set 


ft  »{x6E3:<x-x,Q(x-x)>  -  q13 (x1-x1)2+qJ3x3+qJ3[x3+(oqXi+rqJ3)/(Kq33)]2<a} 

(38) 

which  contains  ft.  Such  a  set  is  given  by  ft*  in  (28)  with 
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(39) 


/a'  -  (|c|/Y)a  +  ||Q2(xg-x)  |j  -  [|cj/(2  min{l,b,o}]^ 

+  (l/2){q11Xj(l  +  ^qjjqjj)  +  (aqu  +  rq,3) 2/ (K2q3 3) }2, 

where  C  is  defined  in  (37).  We  wish  to  choose  qn,q33  and  x’j  such  that 
the  volume  of  ft' 


V  *’  (a')3/qu)2q33  “  3— 3^7q -  j( |c  1  hf*  i[q i £ (: 

(aq1i+rq33)2-]'r)3 
+  — 


1  + 


a2q 


q» 


r) 


K2q33 


rr 


(40) 


is  minimized.  Since  for  any  fixed  q1i»q33>0»  V^,  is  a  strictly  monotone 
increasing  function  of  x3  for  Xj^O,  hence  we  set  Xj-0  in  (40).  The  re¬ 
sulting  V^,  can  be  rewritten  in  terms  of  the  ratio  q*q33/qu 


_ _  _£ 

Vflf  -  (47r/3)t^/(2Y)  +  1/2]  *  | K | — * q  2(o+rq)*, 


(41) 


which  has  a  minimum  point  in  [0,®[  given  by  q  •  5<r/r. 

This  specifies  a  fi'-{xeE1;x?+qx2+q[x3+6r/ (5K)]2<(9/5)  (arK~2)  +  l)2} 

which  contains  the  manifold  generated  by  the  chaotic  oscillations. 

It  is  known  [1]  that  for  K«*-l,  a-10  and  b=8/3,  chaotic  oscillations 
or  a  strange  attractor  exist  for  r > 24.06.  Since  for  K^O,  systems  (2) 
and  (2')  are  equivalent,  we  conclude  that  chaotic  oscillations  exist  for 
the  foregoing  values  of  0,b  and  r  and  any  KfO.  The  peak  magnitudes  | | 
of  the  oscillations  or  the  size  of  the  strange  attractor  are  inversely 
proportional  to  |k|. 

Figure  1  shows  a  typical  buildup  of  the  chaotic  oscillations  in  the 
time  domain  for  the  foregoing  values  of  cr,b,K«l  and  r»30.  The  projections 
of  the  trajectory  onto  the  (x3,x2)  and  (x3,x2)  planes  are  shown  in  Fig. 2. 
For  these  parameter  values,  V^f  given  by  (41)  has  a  minimum  point  q**3/3 
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which  gives  ft*  -  {x£E3:  x2/(61.185)2  +  x2/(47.394)2  +  (x3+36)2/(47.394)2 
^  l}  containing  the  strange  attractor.  The  ellipses  corresponding  to  the 
boundary  of  ft'  in  its  principal-axes  planes  are  also  shown  in  Fig. 2.  It 
should  be  mentioned  that  nonellipsoldal  bounds  for  the  strange  attractor 
of  the  Lorenz  equation  have  been  obtained  using  certain  invariants  of  the 
equation  with  0*0  [13], [14], 

Besides  the  Lorenz  equation,  there  are  a  number  of  other  nonlinear 
systems  which  are  known  to  exhibit  chaotic  oscillations  and  can  be  written 
in  the  form  of  (2).  An  example  is  given  by  Rossler  [15]  where  A,B  and  c 
are 

*  0  -1  -1  To  0  0  ' 

A  =*  1  a  0  ,  B  -  0  0  0  ,  c  -  (1,0, 0)T,  (42) 

.b  0  -rj  L0  0  1. 

and  K— 1.  Here,  A  is  nonsingular  if  0b^r.  The  equilibrium  set  E  con- 

—  1  —1  T 

sists  of  the  origin  and  the  point  (-0b+r,b-rcr  ,ra  -b)  .  It  can  be 

T 

readily  verified  that  (A,B)  is  completely  controllable  and  if  0+ r,  (A,c  ) 
is  completely  observable.  Also, there  does  not  exist  nontrivial  linear 
invariant  subspaces  if  Qj r. 

It  is  known  that  this  system  has  chaotic  oscillations  for  b-0.4, 
r»4.5  and  0.36  ^  a  <  0.5.  Again,  Using  Theorem  2,  we  can 

construct  an  ellipsoidal  set  ft  which  contains  the  chaotic  orbits.  We  omit 
the  details  here.  A  typical  solution  corresponding  to  chaotic  oscillations 
is  shown  in  Figs.  3  and  4. 

VI.  CONCLUDING  REMARKS 

The  existence  of  chaotic  oscillations  in  deterministic  nonlinear 
dynamical  systems  is  an  intriguing  phenomenon  which  is  inherent  in  a  number 
of  mathematical  models  for  real-world  systems.  For  models  of  population 
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dynamics  and  economic  systems ,  chaotic  oscillations  Imply  the  absence  ot 
predictable  "cycles".  In  fluid  models,  such  oscillations  could  provide 
a  mathematical  explanation  of  various  turbulence  phenomena.  Chaotic  oscil¬ 
lations  could  also  arise  in  bilinear  systems  with  linear  feedback  controls 
as  demonstrated  here.  They  cannot  be  predicted  or  analyzed  using  conven¬ 
tional  methods.  Although  at  the  present  time,  there  are  a  number  of  mathe¬ 
matical  results  pertaining  to  various  strange  attractors  in  abstract  dy¬ 
namical  systems  [16], [17].  However,  they  are  not  readily  applicable  for 
determining  the  existence  of  chaotic  oscillations  even  for  the  relatively 
simple  class  of  nonlinear  systems  considered  here.  The  establishment  of 
sufficient  conditions  for  the  existence  of  chaotic  oscillations  or  strange 
attractors  is  a  difficult  challenging  problem  in  nonlinear  system  theory. 

Here,  we  have  considered  only  continuous-time  systems.  Similar  re¬ 
sults  can  be  obtained  for  discrete-time  bilinear  systems  with  linear  feed¬ 
back  controls  of  the  form: 

x(k+l)  -  Ax(k)  +  u(k)Bx(k), 

(43) 

u(k)  -  -K<c,x(k)>. 

It  is  known  that  chaotic  oscillations  can  exist  in  such  systems  also.  In 

2 

fact,  the  scalar  equation:  x(k+l)  =  ax(k)  -  ax  (k)  has  chaotic  oscillations 
for  3.57<a$4.0  [18]. 

A  problem  of  practical  interest  is  the  extension  of  nonturbulent  re¬ 
gime  in  a  fluid  system  to  a  broader  range  of  Reynold's  numbers.  If  fluid 
turbulence  is  indeed  explainable  by  the  theory  of  chaotic  oscillations, 
then  it  may  be  possible  to  achieve  the  desired  result  by  introducing  ap¬ 
propriate  feedback  controls  into  the  system.  Finally,  it  is  of  interest 
to  synthesize  nonlinear  systems  to  generate  chaotic  oscillations  or  pseudo¬ 
random  functions  having  certain  desired  properties.  This  idea  has  al- 
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ready  been  utilized  in  a  number  of  pseudo-random  number  algorithms  in 
digital  computation. 
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FIGURE  CAPTIONS 


Figure  1 

Figure  2 


Figure  3 


Figure  4 


:  Time  domain  solution  of  Lorenz  equation  with  K*l,a»10,b«8/3  and 
r**30;  initial  conditions:  (0)*x2 (0)-50,x3 (0)=0. 

:  Projections  of  the  trajectory  of  Lorenz  equation  (with  parameters 
given  in  Fig.l)  in  the  state  space  onto  the  (x2,x3)  and  (x2,x,) 
planes,  and  the  ellipses  corresponding  co  the  boundary  of  71 '  in 
its  principal-axes  planes  indicated  by  dashed  curves. 

:  Time  domain  solution  of  eq.  (2)  (with  matrices  (42))  for  b=0.*, 
r“4.5,o*4.0  and  K*-l;  initial  conditions:  x3  (0)=x3 (0)»0,  x2(0)» 
4.0. 

:  Projections  of  the  trajectory  of  eq.(2)  (with  matrices  (42)  and 
parameters  given  in  Fig. 3)  in  the  state  space  onto  the  (x2,x3) 
and  (x2,x1)  planes. 
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